rm(list = ls())
Packages included:
# load pacman package from the repository, if you do not already have it
if (!require('pacman')) install.packages('pacman', repos="https://cloud.r-project.org")
pacman::p_load(tidyverse, # set of packages for data manipulation, exploration and visualisation
rjags, # to link JAGS and R
R2jags, # to link JAGS and R
skimr, # for quick dataframe inspection
corrr, # output correlation matrices as data frame
cowplot, # combine plot panels
patchwork, # -"-
rmarkdown) # for R Markdown formatting
Functions used:
# cover plot grids ----
pred.plot.grid <- function(df){
taxa <- df %>% pull(taxon) %>% unique() %>% as.character()
for(taxon_nr in 1:length(taxa)){
plot <- ggplot(data = df %>% filter(taxon == taxa[taxon_nr]),
aes(y = cover, group = taxon)) +
geom_point(aes(x = pred_value), colour = "darkgrey") +
geom_smooth(aes(x = pred_value), method = "lm", colour = "darkgreen", se = TRUE, na.rm = TRUE) +
geom_smooth(aes(x = pred_value), method = "lm", formula = y ~ poly(x, 2), colour = "darkorange", se = TRUE, na.rm = TRUE) +
facet_wrap(~predictor, scales = "free", ncol = 4) +
scale_y_continuous("relative no. pin hits per plot",
limits = c(0, max(df %>%
filter(taxon == taxa[taxon_nr]) %>%
pull(cover)))) +
labs(x = "predictor value") +
ggtitle(paste0(taxa[taxon_nr], " cover ~ predictors")) +
theme_bw() +
theme(legend.position = "none")
print(plot)
}
}
# effect size plots ----
# for cases with only 'significant' effects
model_plot_sig_function <- function(model_coeff_output, title_string, plot_width) {
target_vars <- c("b.tempjja.x", "b.tempjja.x2",
"b.tempcont.x", "b.tempcont.x2",
"b.precipjja.x", "b.precipjja.x2",
"b.sri",
"b.tri",
"b.twi",
"b.compet")
solutions <- model_coeff_output
names(solutions) <- c("variable", "post.mean", "post.sd", "l95", "l90", "u90", "u95", "Rhat")
solutions <- solutions %>%
filter(variable %in% target_vars)
# solutions$variable <- factor(solutions$variable,
# levels = c("b.tempjja.x", "b.tempjja.x2",
# "b.tempcont.x", "b.tempcont.x2",
# "b.precipjja.x", "b.precipjja.x2",
# "b.sri",
# "b.tri",
# "b.twi",
# "b.compet"))
min_value <- floor(min(solutions$l95))
max_value <- ceiling(max(solutions$u95))
solutions$sig <- "ns"
solutions$sig[solutions$l95 < 0 & solutions$u95 < 0] <- "sig"
solutions$sig[solutions$l95 > 0 & solutions$u95 > 0] <- "sig"
label_colour <- rep("black", nrow(solutions))
label_colour[solutions$sig == "sig"] <- theme_darkgreen
label_face <- rep("plain", nrow(solutions))
label_face[solutions$sig == "sig"] <- "bold"
# label_face[response == "T1_mean" & solutions$sig == "sig"] <- "bold"
title_string <- title_string
title_colour <- "grey10"
# if(response == "T1_mean" | response == "T1_amp") title_colour <- theme_red
# if(response == "T2_mean" | response == "T2_amp") title_colour <- theme_yellow
# if(response == "T1_mean") response <- "Soil"
# if(response == "T2_mean") response <- "Ground"
model_plot_sig <- ggplot(solutions, aes(x = variable, y = post.mean,
ymin = l95, ymax = u95,
colour = sig)) +
geom_point() +
geom_errorbar(width = .8) +
theme_cowplot(18) +
ylab("Effect Size (scaled)") +
xlab("") +
ggtitle(paste0(title_string)) +
scale_colour_manual(values = c("black", theme_darkgreen)) +
scale_y_continuous(limits = c(min_value, max_value), breaks = seq(min_value,max_value,0.5)) +
# scale_x_discrete(limits = c("b.tempjja.x", "b.tempjja.x2",
# "b.tempcont.x", "b.tempcont.x2",
# "b.precipjja.x", "b.precipjja.x2",
# "b.sri",
# "b.tri",
# "b.twi",
# "b.compet"),
# labels = c("summer temperature", bquote(.("summer") *" "* temperature^2),
# "temperature variability", bquote(.("temperature") *" "* variability^2),
# "summer precipitation", bquote(.("summer") *" "* precipitation^2),
# "solar radiation",
# "terrain ruggedness",
# "moisture availability",
# "competition")) +
annotate("segment", x = 0, xend = plot_width, y = 0, yend = 0) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = label_colour, face = label_face),
plot.title = element_text(colour = title_colour, face = "italic"),
legend.position = "none")
return(model_plot_sig)
}
# for cases with marginal 'significance'
model_plot_marg_function <- function(model_coeff_output, title_string, plot_width) {
target_vars <- c("b.tempjja.x", "b.tempjja.x2",
"b.tempcont.x", "b.tempcont.x2",
"b.precipjja.x", "b.precipjja.x2",
"b.sri",
"b.tri",
"b.twi",
"b.compet")
solutions <- model_coeff_output
names(solutions) <- c("variable", "post.mean", "post.sd", "l95", "l90", "u90", "u95", "Rhat")
solutions <- solutions %>%
filter(variable %in% target_vars)
# solutions$variable <- factor(solutions$variable,
# levels = c("b.tempjja.x", "b.tempjja.x2",
# "b.tempcont.x", "b.tempcont.x2",
# "b.precipjja.x", "b.precipjja.x2",
# "b.sri",
# "b.tri",
# "b.twi",
# "b.compet"))
min_value <- floor(min(solutions$l95))
max_value <- ceiling(max(solutions$u95))
solutions$sig <- "ns"
solutions$sig[solutions$l95 < 0 & solutions$u95 < 0] <- "sig"
solutions$sig[solutions$l95 > 0 & solutions$u95 > 0] <- "sig"
solutions$sig[solutions$l90 < 0 & solutions$u90 < 0 & solutions$l95 < 0 & solutions$u95 > 0] <- "marg"
solutions$sig[solutions$l90 > 0 & solutions$u90 > 0 & solutions$l95 < 0 & solutions$u95 > 0] <- "marg"
label_colour <- rep("black", nrow(solutions))
label_colour[solutions$sig == "sig"] <- theme_darkgreen
label_colour[solutions$sig == "marg"] <- theme_purple
label_face <- rep("plain", nrow(solutions))
label_face[solutions$sig == "sig"] <- "bold"
# label_face[response == "T1_mean" & solutions$sig == "sig"] <- "bold"
title_string <- title_string
title_colour <- "grey10"
# if(response == "T1_mean" | response == "T1_amp") title_colour <- theme_red
# if(response == "T2_mean" | response == "T2_amp") title_colour <- theme_yellow
# if(response == "T1_mean") response <- "Soil"
# if(response == "T2_mean") response <- "Ground"
model_plot_marg <- ggplot(solutions, aes(x = variable, y = post.mean,
ymin = l95, ymax = u95,
colour = sig)) +
geom_point() +
geom_errorbar(width = .8) +
theme_cowplot(18) +
ylab("Effect Size (scaled)") +
xlab("") +
ggtitle(paste0(title_string)) +
scale_colour_manual(values = c(theme_purple, "black", theme_darkgreen)) +
scale_y_continuous(limits = c(min_value, max_value), breaks = seq(min_value,max_value,0.5)) +
# scale_x_discrete(limits = c("b.tempjja.x", "b.tempjja.x2",
# "b.tempcont.x", "b.tempcont.x2",
# "b.precipjja.x", "b.precipjja.x2",
# "b.sri",
# "b.tri",
# "b.twi",
# "b.compet"),
# labels = c("summer temperature", bquote(.("summer") *" "* temperature^2),
# "temperature variability", bquote(.("temperature") *" "* variability^2),
# "summer precipitation", bquote(.("summer") *" "* precipitation^2),
# "solar radiation",
# "terrain ruggedness",
# "moisture availability",
# "competition")) +
annotate("segment", x = 0, xend = plot_width, y = 0, yend = 0) +
theme(axis.text.x = element_text(angle = 45, hjust = 1, colour = label_colour, face = label_face),
plot.title = element_text(colour = title_colour, face = "italic"),
legend.position = "none")
return(model_plot_marg)
}
Colour scheme:
# Colour scheme ####
theme_red <- "#E64447FF"
theme_blue <- "#12435DFF"
theme_lightblue <- "#1FA3AEFF"
theme_grey <- "#EDEDEDFF"
theme_yellow <- "#EEBE5BFF"
# additional colours
theme_darkblue <- "#1D5799"
theme_darkgreen <- "#13944D"
theme_orange <- "#B56A24"
theme_purple <- "#8757B3"
This dataset is a new one containing predictors extracted from downscaled CHELSA climate data.
Analyses conducted on fusion table at plot \(\times\) taxon level (as we are explicitly interested in plot-level variation of predictors, esp. slope, solar radiation, terrain ruggedness, topographical wetness, and competition): 3726 observations for 38 variables
env_cov_bio <- read.csv("../data/nuuk_env_cover_plots.csv", header = T) %>%
# rename variable "twi_90m" to "twi"
rename(twi = twi_90m)
including predictors
and response variable
env_cov_bio_sub <- env_cov_bio %>%
select(site_alt_plotgroup_id, plot, site, site_alt_id, year, long, lat, alt, # plot info / metadata
ends_with("_ts_30"), # CHELSA predictors averaged over 10-year period prior to study year
inclin_down, twi, tri, sri,
#mean_summ_ndvi_yos, cv_mean_summ_ndvi_2001_to_yos, Perc_dist_coast_lines, # environmental data
taxon, cover, compet) # taxon, cover response, competition pressure
head(env_cov_bio_sub)
Let’s check for correlation between the different moisture predictors and terrain ruggedness:
(cor_moist <- env_cov_bio %>%
select(twi, ndwi, tcws, tri)) %>%
correlate(diagonal = 1)
Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'
\(\Rightarrow\) NDWI and TCWS are highly correlated, but both are largely independent from TWI and TRI. We’ll go with TWI for now, as it is least confounded by vegetation, but will have to discuss maybe including a second index, depending on hypothesis/ecological implication.
Predictors don’t always vary between plots within plotgroups - perhaps due to several falling into the same CHELSA grid cell.
Example: plots within site 1, altitude 20, plot group 1: plot P146 is slightly off and therefore has different climate variables than the other ones
env_cov_bio %>% filter(taxon == "Betula nana" & site_alt_plotgroup_id == "1_20_1") %>%
ggplot(aes(x = long, y = lat)) +
geom_point() +
geom_text(aes(label = plot), hjust = 0.0001) +
xlim(c(-51.78675, -51.7863))
env_cov_bio %>% filter(taxon == "Betula nana" & site_alt_plotgroup_id == "1_20_1") %>%
select(site_alt_plotgroup_id, plot, tempjja_ts_30)
Check for correlation between predictors:
predictors_set <- env_cov_bio_sub %>%
select(ends_with("_ts_30"), # CHELSA predictors averaged over 10-year period prior to study year
inclin_down, sri, tri, twi, # environmental data
compet) %>%
names()
# create basic correlation matrix
(cor_mat <- env_cov_bio_sub %>%
dplyr::select(predictors_set) %>%
correlate(diagonal = 1) %>%
# drop all values < .4 to increase readability
mutate_if(is.numeric, ~ ifelse(abs(.) < .4, NA, .)))
Correlation method: 'pearson'
Missing treated using: 'pairwise.complete.obs'
\(\Rightarrow\) Maximum temperature is highly correlated with JJA temperature, minimum temperature with continentality, slope (inclin_down) with solar radiation, and spring precipitation variables with JJA precipitation. Let’s exclude them step by step and check the variance inflation factors (VIF) along the way. VIF values > 5 indicate collinearity issues.
# for the whole set of predictors
(vif_predictors_1 <- env_cov_bio_sub %>%
dplyr::select(predictors_set) %>%
usdm::vif()) #%>% View()
# => drop tempmax (correlated w/ tempjja)
(vif_predictors_2 <- env_cov_bio_sub %>%
dplyr::select(predictors_set,
-tempmax_ts_30) %>%
usdm::vif()) #%>% View()
# => drop tempmin (correlated w/ tempcont)
(vif_predictors_3 <- env_cov_bio_sub %>%
dplyr::select(predictors_set,
-tempmax_ts_30,
-tempmin_ts_30) %>%
usdm::vif()) #%>% View()
# => drop precipjfmam (correlated w/ precipjja & tempcont)
(vif_predictors_4 <- env_cov_bio_sub %>%
dplyr::select(predictors_set,
-tempmax_ts_30,
-tempmin_ts_30,
-contains("jfm")) %>%
usdm::vif()) #%>% View()
# => drop inclin_down (correlated w/ SRI)
(vif_predictors_5 <- env_cov_bio_sub %>%
dplyr::select(predictors_set,
-tempmax_ts_30,
-tempmin_ts_30,
-inclin_down,
-contains("jfm")) %>%
usdm::vif()) #%>% View()
(vif_predictors_6 <- env_cov_bio_sub %>%
dplyr::select(predictors_set,
-tempmax_ts_30,
-tempmin_ts_30,
-inclin_down,
-contains("mam")) %>%
usdm::vif()) #%>% View()
All VIF are now < 3.1 -> OK! We can now exclude the dropped variables to obtain the final dataset:
env_cov_bio_sub <- env_cov_bio_sub %>%
dplyr::select(-tempmax_ts_30,
-tempmin_ts_30,
-inclin_down,
-contains("mam"))
(nuuk_spec_abundance_plot <- ggplot(env_cov_bio %>%
# make site a factor
mutate(site = as.factor(site)) %>%
mutate(site_alt_id = factor(site_alt_id, levels = c(paste(rep(1, 3), c("20", "100", "200"), sep = "_"),
paste(rep(2, 3), c("20", "100", "200"), sep = "_"),
paste(rep(3, 5), c("20", "100", "200", "300", "400"), sep = "_"),
paste(rep(4, 6), c("20", "100", "200", "300", "400", "500"), sep = "_"),
paste(rep(5, 6), c("20", "100", "200", "300", "400", "500"), sep = "_")))) %>%
# group by site and isocline
group_by(site, site_alt_id),
aes(x = site_alt_id,
y = cover,
fill = site)) +
# draw boxplots of cover
geom_boxplot() +
# split by taxon
facet_grid(rows = vars(taxon)) +
# scale_fill_manual() +
theme_bw() +
xlab("site / isocline") +
theme(strip.text.y = element_text(face = "italic"),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = rel(1.2)),
axis.title = element_text(size = rel(1.2))))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_species_abundance_gradient.eps"),
# nuuk_spec_abundance_plot, base_height = 18, base_aspect_ratio = 0.8)
# data compilation ----
predictors_set <- env_cov_bio %>%
select(ends_with("_ts_30"), # CHELSA predictors averaged over 10-year period prior to study year
inclin_down, sri, tri, twi, # environmental data
compet) %>%
names()
preds_plot_data <- env_cov_bio %>%
# convert site, isocline & plotgroup to factors
mutate(site = as.factor(site)) %>%
mutate(site_alt_id = factor(site_alt_id, levels = c(paste(rep(1, 3), c("20", "100", "200"), sep = "_"),
paste(rep(2, 3), c("20", "100", "200"), sep = "_"),
paste(rep(3, 5), c("20", "100", "200", "300", "400"), sep = "_"),
paste(rep(4, 6), c("20", "100", "200", "300", "400", "500"), sep = "_"),
paste(rep(5, 6), c("20", "100", "200", "300", "400", "500"), sep = "_")))) %>%
mutate(plotgroup = as.factor(plotgroup)) %>%
select(site, site_alt_id, predictors_set) %>%
# pivot to long format
pivot_longer(cols = predictors_set,
names_to = "predictor",
values_to = "value") %>%
# convert predictor col to factor
mutate(predictor = as.factor(predictor)) %>%
# rename predictors
mutate(predictor = fct_recode(predictor,
"summer \n temperature [°C]" = "tempjja_ts_30",
"yearly maximum \n temperature [°C]" = "tempmax_ts_30",
"yearly minimum \n temperature [°C]" = "tempmin_ts_30",
"annual temperature \n variability [°C]" = "tempcont_ts_30",
"cumulative summer \n precipitation [mm]" = "precipjja_ts_30",
"cumulative winter-spring \n precipitation [mm]" = "precipjfmam_ts_30",
"cumulative spring \n precipitation [mm]" = "precipmam_ts_30",
"slope angle [°]" = "inclin_down",
"Solar Radiation \n Index" = "sri",
"Terrain Ruggedness \n Index" = "tri",
"Topographic Wetness \n Index" = "twi",
"overgrowing \n competition" = "compet"),
predictor = fct_relevel(predictor,
"summer \n temperature [°C]",
"yearly maximum \n temperature [°C]",
"yearly minimum \n temperature [°C]",
"annual temperature \n variability [°C]",
"cumulative summer \n precipitation [mm]",
"cumulative winter-spring \n precipitation [mm]",
"cumulative spring \n precipitation [mm]",
"slope angle [°]",
"Solar Radiation \n Index",
"Terrain Ruggedness \n Index",
"Topographic Wetness \n Index",
"overgrowing \n competition")) %>%
# group by site & isocline
group_by(site, site_alt_id)
# climatic predictors ----
predictors_set_clim_long <- c("summer \n temperature [°C]",
"yearly maximum \n temperature [°C]",
"yearly minimum \n temperature [°C]",
"annual temperature \n variability [°C]",
"cumulative summer \n precipitation [mm]",
"cumulative winter-spring \n precipitation [mm]",
"cumulative spring \n precipitation [mm]")
(nuuk_preds_clim_gradient_plot <- ggplot(preds_plot_data %>%
filter(predictor %in% predictors_set_clim_long),
aes(x = site_alt_id,
y = value,
fill = site)) +
# draw boxplots of values
geom_boxplot() +
# split by taxon
facet_grid(rows = vars(predictor), scales = "free_y") +
# scale_fill_manual() +
theme_bw() +
xlab("site / isocline") +
theme(strip.text.y = element_text(face = "italic"),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = rel(1.2)),
axis.title = element_text(size = rel(1.2))))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_preds_clim_gradient.eps"),
# nuuk_preds_clim_gradient_plot, base_height = 15, base_aspect_ratio = 0.8)
# environmental predictors ----
predictors_set_env_long <- c("slope angle [°]",
"Solar Radiation \n Index",
"Terrain Ruggedness \n Index",
"Topographic Wetness \n Index",
"overgrowing \n competition")
(nuuk_preds_env_gradient_plot <- ggplot(preds_plot_data %>%
filter(predictor %in% predictors_set_env_long),
aes(x = site_alt_id,
y = value,
fill = site)) +
# draw boxplots of values
geom_boxplot() +
# split by taxon
facet_grid(rows = vars(predictor), scales = "free_y") +
# scale_fill_manual() +
theme_bw() +
xlab("site / isocline") +
theme(strip.text.y = element_text(face = "italic"),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = rel(1.2)),
axis.title = element_text(size = rel(1.2))))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_preds_env_gradient.eps"),
# nuuk_preds_env_gradient_plot, base_height = 15, base_aspect_ratio = 0.8)
# final set of predictors ----
predictors_set_final_long <- c("summer \n temperature [°C]",
"annual temperature \n variability [°C]",
"cumulative summer \n precipitation [mm]",
"Solar Radiation \n Index",
"Terrain Ruggedness \n Index",
"Topographic Wetness \n Index",
"overgrowing \n competition")
(nuuk_preds_final_gradient_plot <- ggplot(preds_plot_data %>%
filter(predictor %in% predictors_set_final_long),
aes(x = site_alt_id,
y = value,
fill = site)) +
# draw boxplots of values
geom_boxplot() +
# split by taxon
facet_grid(rows = vars(predictor), scales = "free_y") +
# scale_fill_manual() +
theme_bw() +
xlab("site / isocline") +
theme(strip.text.y = element_text(face = "italic"),
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1, size = rel(1.2)),
axis.title = element_text(size = rel(1.2))))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_preds_final_gradient.eps"),
# nuuk_preds_final_gradient_plot, base_height = 15, base_aspect_ratio = 0.8)
How do the trends in cover response look for each of the predictors? Plot panels below show linear (green) and quadratic (orange) trendlines derived from linear models (different than the Bayesian approach used for analysis! but anyway useful to catch trends):
# Prepare data for panel plots for cover against predictors for each species
env_cov_bio.long <- env_cov_bio %>%
select(taxon,
tempjja_ts_30,
tempcont_ts_30,
precipjja_ts_30,
sri,
tri,
twi,
compet,
cover) %>%
# pivot to long format
pivot_longer(cols = c(tempjja_ts_30,
tempcont_ts_30,
precipjja_ts_30,
sri,
tri,
twi,
compet),
names_to = "predictor",
values_to = "pred_value")
# run
pred.plot.grid(env_cov_bio.long)
Data was ordered by site/altitude/plotgroup and taxon
env_cov_bio_sub <- env_cov_bio_sub[order(env_cov_bio_sub$site_alt_plotgroup_id, env_cov_bio_sub$taxon),]
As JAGS is only able to handle numeric input, all variables are assigned a numeric identifier:
env_cov_bio_sub$plotgroup.NUM <- as.numeric(factor(env_cov_bio_sub$site_alt_plotgroup_id,
levels = unique(env_cov_bio_sub$site_alt_plotgroup_id)))
env_cov_bio_sub$plot.NUM <- as.numeric(factor(env_cov_bio_sub$plot,
levels = unique(env_cov_bio_sub$plot)))
env_cov_bio_sub$site_alt.NUM <- as.numeric(factor(env_cov_bio_sub$site_alt_id,
levels = unique(env_cov_bio_sub$site_alt_id)))
env_cov_bio_sub$site.NUM <- as.numeric(factor(env_cov_bio_sub$site, levels = unique(env_cov_bio_sub$site)))
env_cov_bio_sub$taxon.NUM <- as.numeric(factor(env_cov_bio_sub$taxon, levels = unique(env_cov_bio_sub$taxon)))
Taxa were coded as follows:
data.frame(taxon = levels(env_cov_bio_sub$taxon),
num = unique(env_cov_bio_sub$taxon.NUM))
Numeric predictors were scaled and centered:
num_pred <- env_cov_bio_sub %>% select(alt,
ends_with("_ts_30"),
sri,
tri,
starts_with("twi"),
matches("compet"))
for(i in 1:length(num_pred)){
col <- colnames(num_pred[i])
env_cov_bio_sub[paste0(col,"C")] <- as.numeric(scale(num_pred[i], scale = TRUE, center = TRUE))
}
To account for the range of the cover response (\(0 \leq cover \leq 1\)), the model needs a mixed structure incorporating a beta distribution (for all continuous values with \(0 < cover < 1\)) and a binomial distribution (for all discrete values of \(cover = \{0, 1\}\)). An additional variable cover_discrete was introduced to separate the dataset into discrete (= 1) and continuous (= 0) cover values:
env_cov_bio_sub$cover_discrete <- ifelse(env_cov_bio_sub$cover == 1 | env_cov_bio_sub$cover == 0, 1, 0)
The dataset was then ready to be split up into the species of interest. As a first trial, I focused on Betula nana:
JAGS needs data input in list format, so I provided all relevant variables as follows:
…and the parameters to be monitored:
Notes at first glance:
As the model converges well, we can proceed to apply it to all of the species.
The dataset was then ready to be split up into the species of interest. We create separate data subsets for all/discrete/continuous response variable values for each species:
# split dataframe by taxon
env_cov_bio_sub_spec.tot <- split(env_cov_bio_sub, env_cov_bio_sub$taxon)
# assign taxon name to list elements
# >> for total datasets
for (taxon_id in 1:nlevels(env_cov_bio_sub$taxon)){
# extract 3-letter genus name string
assign(paste0(str_extract(levels(env_cov_bio_sub$taxon)[taxon_id],
"^\\w{3}"),
# extract and capitalise 3-letter species name string
str_to_title(str_remove(str_extract(levels(env_cov_bio_sub$taxon)[taxon_id],
"\\s\\w{3}"),
"\\s")),
# add extension
".tot"),
# assign to respective list element
env_cov_bio_sub_spec.tot[[taxon_id]])
}
# >> for discrete datasets
env_cov_bio_sub_spec.dis <- list()
for (taxon_id in 1:nlevels(env_cov_bio_sub$taxon)){
# filter for discrete response values
env_cov_bio_sub_spec.dis[[taxon_id]] <- filter(env_cov_bio_sub_spec.tot[[taxon_id]], cover_discrete == 1)
# extract 3-letter genus name string
assign(paste0(str_extract(levels(env_cov_bio_sub$taxon)[taxon_id],
"^\\w{3}"),
# extract and capitalise 3-letter species name string
str_to_title(str_remove(str_extract(levels(env_cov_bio_sub$taxon)[taxon_id],
"\\s\\w{3}"),
"\\s")),
# add extension
".dis"),
# assign to respective list element
env_cov_bio_sub_spec.dis[[taxon_id]])
}
# >> for continuous datasets
env_cov_bio_sub_spec.cont <- list()
for (taxon_id in 1:nlevels(env_cov_bio_sub$taxon)){
# filter for continuous response values
env_cov_bio_sub_spec.cont[[taxon_id]] <- filter(env_cov_bio_sub_spec.tot[[taxon_id]], cover_discrete == 0)
# extract 3-letter genus name string
assign(paste0(str_extract(levels(env_cov_bio_sub$taxon)[taxon_id],
"^\\w{3}"),
# extract and capitalise 3-letter species name string
str_to_title(str_remove(str_extract(levels(env_cov_bio_sub$taxon)[taxon_id],
"\\s\\w{3}"),
"\\s")),
# add extension
".cont"),
# assign to respective list element
env_cov_bio_sub_spec.cont[[taxon_id]])
}
JAGS needs data input in list format, so I provided all relevant variables as follows:
# Compile data into lists
# BetNan ----
shrub_gradient_jags.BetNan.data <- list(
# plot level predictors, for discrete...
cov.dis = BetNan.dis$cover,
plotgroup.dis = BetNan.dis$plotgroup.NUM, #AB added this
# isocline.dis = BetNan.dis$site_alt.NUM,
# inclin_down.dis = BetNan.dis$inclin_downC,
sri.dis = BetNan.dis$sriC,
tri.dis = BetNan.dis$triC,
twi.dis = BetNan.dis$twiC,
compet.dis = BetNan.dis$competC,
N_discrete = nrow(BetNan.dis),
# ...and continuous part of the data
cov.cont = BetNan.cont$cover,
plotgroup.cont = BetNan.cont$plotgroup.NUM, #AB added this
# isocline.cont = BetNan.cont$site_alt.NUM,
# inclin_down.cont = BetNan.cont$inclin_downC,
sri.cont = BetNan.cont$sriC,
tri.cont = BetNan.cont$triC,
twi.cont = BetNan.cont$twiC,
compet.cont = BetNan.cont$competC,
N_cont = nrow(BetNan.cont),
# plot group level predictors
tempjja.tot = BetNan.tot$tempjja_ts_30C[!duplicated(BetNan.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = BetNan.tot$tempmax_ts_30C[!duplicated(BetNan.tot$plotgroup.NUM)],
# tempmin.tot = BetNan.tot$tempmin_ts_30C[!duplicated(BetNan.tot$plotgroup.NUM)],
tempcont.tot = BetNan.tot$tempcont_ts_30C[!duplicated(BetNan.tot$plotgroup.NUM)],
precipjja.tot = BetNan.tot$precipjja_ts_30C[!duplicated(BetNan.tot$plotgroup.NUM)],
# precipjfmam.tot = BetNan.tot$precipjfmam_ts_30C[!duplicated(BetNan.tot$plotgroup.NUM)]
N_plotgroups = length(unique(BetNan.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = BetNan.tot$altC[!duplicated(BetNan.tot$site_alt.NUM)],
# N_isoclines = length(unique(BetNan.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(BetNan.tot$competC), to = max(BetNan.tot$competC), length.out = 100),
xhat_sri = seq(from = min(BetNan.tot$sriC), to = max(BetNan.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(BetNan.tot$triC), to = max(BetNan.tot$triC), length.out = 100),
xhat_twi = seq(from = min(BetNan.tot$twiC), to = max(BetNan.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(BetNan.tot$tempjja_ts_30C), to = max(BetNan.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(BetNan.tot$precipjja_ts_30C), to = max(BetNan.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(BetNan.tot$tempcont_ts_30C), to = max(BetNan.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(BetNan.tot$tempjja_ts_30C,0.05),quantile(BetNan.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.BetNan.data)
List of 28
$ cov.dis : num [1:223] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:223] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.dis : num [1:223] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.dis : num [1:223] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.dis : num [1:223] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.dis : num [1:223] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_discrete : int 223
$ cov.cont : num [1:191] 0.08 0.04 0.08 0.12 0.08 0.56 0.24 0.2 0.08 0.12 ...
$ plotgroup.cont: num [1:191] 5 5 6 10 10 11 11 11 11 12 ...
$ sri.cont : num [1:191] 0.823 -0.248 -0.663 -0.678 -0.192 ...
$ tri.cont : num [1:191] 0.4613 0.3752 -0.0373 -0.2759 0.1133 ...
$ twi.cont : num [1:191] 1.78 1.78 1.22 -1.29 -1.29 ...
$ compet.cont : num [1:191] -0.477 -0.708 -0.708 -0.708 -0.708 ...
$ N_cont : int 191
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.671 -0.633 -0.596 -0.559 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# CasTet ----
shrub_gradient_jags.CasTet.data <- list(
# plot level predictors, for discrete...
cov.dis = CasTet.dis$cover,
plotgroup.dis = CasTet.dis$plotgroup.NUM, #AB added this
# isocline.dis = CasTet.dis$site_alt.NUM,
# inclin_down.dis = CasTet.dis$inclin_downC,
sri.dis = CasTet.dis$sriC,
tri.dis = CasTet.dis$triC,
twi.dis = CasTet.dis$twiC,
compet.dis = CasTet.dis$competC,
N_discrete = nrow(CasTet.dis),
# ...and continuous part of the data
cov.cont = CasTet.cont$cover,
plotgroup.cont = CasTet.cont$plotgroup.NUM, #AB added this
# isocline.cont = CasTet.cont$site_alt.NUM,
# inclin_down.cont = CasTet.cont$inclin_downC,
sri.cont = CasTet.cont$sriC,
tri.cont = CasTet.cont$triC,
twi.cont = CasTet.cont$twiC,
compet.cont = CasTet.cont$competC,
N_cont = nrow(CasTet.cont),
# plot group level predictors
tempjja.tot = CasTet.tot$tempjja_ts_30C[!duplicated(CasTet.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = CasTet.tot$tempmax_ts_30C[!duplicated(CasTet.tot$plotgroup.NUM)],
# tempmin.tot = CasTet.tot$tempmin_ts_30C[!duplicated(CasTet.tot$plotgroup.NUM)],
tempcont.tot = CasTet.tot$tempcont_ts_30C[!duplicated(CasTet.tot$plotgroup.NUM)],
precipjja.tot = CasTet.tot$precipjja_ts_30C[!duplicated(CasTet.tot$plotgroup.NUM)],
# precipjfmam.tot = CasTet.tot$precipjfmam_ts_30C[!duplicated(CasTet.tot$plotgroup.NUM)]
N_plotgroups = length(unique(CasTet.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = CasTet.tot$altC[!duplicated(CasTet.tot$site_alt.NUM)],
# N_isoclines = length(unique(CasTet.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(CasTet.tot$competC), to = max(CasTet.tot$competC), length.out = 100),
xhat_sri = seq(from = min(CasTet.tot$sriC), to = max(CasTet.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(CasTet.tot$triC), to = max(CasTet.tot$triC), length.out = 100),
xhat_twi = seq(from = min(CasTet.tot$twiC), to = max(CasTet.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(CasTet.tot$tempjja_ts_30C), to = max(CasTet.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(CasTet.tot$precipjja_ts_30C), to = max(CasTet.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(CasTet.tot$tempcont_ts_30C), to = max(CasTet.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(CasTet.tot$tempjja_ts_30C,0.05),quantile(CasTet.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.CasTet.data)
List of 28
$ cov.dis : num [1:398] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:398] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.dis : num [1:398] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.dis : num [1:398] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.dis : num [1:398] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.dis : num [1:398] 0.796 1.722 1.028 1.375 0.796 ...
$ N_discrete : int 398
$ cov.cont : num [1:16] 0.04 0.04 0.08 0.04 0.08 0.72 0.08 0.04 0.16 0.04 ...
$ plotgroup.cont: num [1:16] 30 30 30 38 43 49 49 49 49 49 ...
$ sri.cont : num [1:16] -0.314 -1.308 -0.87 0.754 -0.746 ...
$ tri.cont : num [1:16] 1.769 1.565 1.324 -0.807 1.567 ...
$ twi.cont : num [1:16] -1.36 -1.36 -1.31 1.9 -1.2 ...
$ compet.cont : num [1:16] 2.3011 1.1437 2.4168 -0.0138 0.5649 ...
$ N_cont : int 16
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.656 -0.603 -0.55 -0.498 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# EmpNig ----
shrub_gradient_jags.EmpNig.data <- list(
# plot level predictors, for discrete...
cov.dis = EmpNig.dis$cover,
plotgroup.dis = EmpNig.dis$plotgroup.NUM, #AB added this
# isocline.dis = EmpNig.dis$site_alt.NUM,
# inclin_down.dis = EmpNig.dis$inclin_downC,
sri.dis = EmpNig.dis$sriC,
tri.dis = EmpNig.dis$triC,
twi.dis = EmpNig.dis$twiC,
compet.dis = EmpNig.dis$competC,
N_discrete = nrow(EmpNig.dis),
# ...and continuous part of the data
cov.cont = EmpNig.cont$cover,
plotgroup.cont = EmpNig.cont$plotgroup.NUM, #AB added this
# isocline.cont = EmpNig.cont$site_alt.NUM,
# inclin_down.cont = EmpNig.cont$inclin_downC,
sri.cont = EmpNig.cont$sriC,
tri.cont = EmpNig.cont$triC,
twi.cont = EmpNig.cont$twiC,
compet.cont = EmpNig.cont$competC,
N_cont = nrow(EmpNig.cont),
# plot group level predictors
tempjja.tot = EmpNig.tot$tempjja_ts_30C[!duplicated(EmpNig.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = EmpNig.tot$tempmax_ts_30C[!duplicated(EmpNig.tot$plotgroup.NUM)],
# tempmin.tot = EmpNig.tot$tempmin_ts_30C[!duplicated(EmpNig.tot$plotgroup.NUM)],
tempcont.tot = EmpNig.tot$tempcont_ts_30C[!duplicated(EmpNig.tot$plotgroup.NUM)],
precipjja.tot = EmpNig.tot$precipjja_ts_30C[!duplicated(EmpNig.tot$plotgroup.NUM)],
# precipjfmam.tot = EmpNig.tot$precipjfmam_ts_30C[!duplicated(EmpNig.tot$plotgroup.NUM)]
N_plotgroups = length(unique(EmpNig.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = EmpNig.tot$altC[!duplicated(EmpNig.tot$site_alt.NUM)],
# N_isoclines = length(unique(EmpNig.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(EmpNig.tot$competC), to = max(EmpNig.tot$competC), length.out = 100),
xhat_sri = seq(from = min(EmpNig.tot$sriC), to = max(EmpNig.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(EmpNig.tot$triC), to = max(EmpNig.tot$triC), length.out = 100),
xhat_twi = seq(from = min(EmpNig.tot$twiC), to = max(EmpNig.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(EmpNig.tot$tempjja_ts_30C), to = max(EmpNig.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(EmpNig.tot$precipjja_ts_30C), to = max(EmpNig.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(EmpNig.tot$tempcont_ts_30C), to = max(EmpNig.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(EmpNig.tot$tempjja_ts_30C,0.05),quantile(EmpNig.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.EmpNig.data)
List of 28
$ cov.dis : num [1:229] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:229] 2 2 3 5 5 6 6 6 7 7 ...
$ sri.dis : num [1:229] 0.9453 -1.4347 -1.4755 0.7836 -0.0244 ...
$ tri.dis : num [1:229] -0.9834 0.1575 0.0119 0.1007 0.1954 ...
$ twi.dis : num [1:229] -0.754 -0.73 -1.069 1.775 2.121 ...
$ compet.dis : num [1:229] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_discrete : int 229
$ cov.cont : num [1:185] 0.52 0.84 0.6 0.64 0.52 0.16 0.16 0.64 0.16 0.52 ...
$ plotgroup.cont: num [1:185] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.cont : num [1:185] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.cont : num [1:185] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.cont : num [1:185] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.cont : num [1:185] -0.708 -0.708 -0.708 -0.477 -0.708 ...
$ N_cont : int 185
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.665 -0.622 -0.578 -0.535 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# PhyCae ----
shrub_gradient_jags.PhyCae.data <- list(
# plot level predictors, for discrete...
cov.dis = PhyCae.dis$cover,
plotgroup.dis = PhyCae.dis$plotgroup.NUM, #AB added this
# isocline.dis = PhyCae.dis$site_alt.NUM,
# inclin_down.dis = PhyCae.dis$inclin_downC,
sri.dis = PhyCae.dis$sriC,
tri.dis = PhyCae.dis$triC,
twi.dis = PhyCae.dis$twiC,
compet.dis = PhyCae.dis$competC,
N_discrete = nrow(PhyCae.dis),
# ...and continuous part of the data
cov.cont = PhyCae.cont$cover,
plotgroup.cont = PhyCae.cont$plotgroup.NUM, #AB added this
# isocline.cont = PhyCae.cont$site_alt.NUM,
# inclin_down.cont = PhyCae.cont$inclin_downC,
sri.cont = PhyCae.cont$sriC,
tri.cont = PhyCae.cont$triC,
twi.cont = PhyCae.cont$twiC,
compet.cont = PhyCae.cont$competC,
N_cont = nrow(PhyCae.cont),
# plot group level predictors
tempjja.tot = PhyCae.tot$tempjja_ts_30C[!duplicated(PhyCae.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = PhyCae.tot$tempmax_ts_30C[!duplicated(PhyCae.tot$plotgroup.NUM)],
# tempmin.tot = PhyCae.tot$tempmin_ts_30C[!duplicated(PhyCae.tot$plotgroup.NUM)],
tempcont.tot = PhyCae.tot$tempcont_ts_30C[!duplicated(PhyCae.tot$plotgroup.NUM)],
precipjja.tot = PhyCae.tot$precipjja_ts_30C[!duplicated(PhyCae.tot$plotgroup.NUM)],
# precipjfmam.tot = PhyCae.tot$precipjfmam_ts_30C[!duplicated(PhyCae.tot$plotgroup.NUM)]
N_plotgroups = length(unique(PhyCae.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = PhyCae.tot$altC[!duplicated(PhyCae.tot$site_alt.NUM)],
# N_isoclines = length(unique(PhyCae.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(PhyCae.tot$competC), to = max(PhyCae.tot$competC), length.out = 100),
xhat_sri = seq(from = min(PhyCae.tot$sriC), to = max(PhyCae.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(PhyCae.tot$triC), to = max(PhyCae.tot$triC), length.out = 100),
xhat_twi = seq(from = min(PhyCae.tot$twiC), to = max(PhyCae.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(PhyCae.tot$tempjja_ts_30C), to = max(PhyCae.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(PhyCae.tot$precipjja_ts_30C), to = max(PhyCae.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(PhyCae.tot$tempcont_ts_30C), to = max(PhyCae.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(PhyCae.tot$tempjja_ts_30C,0.05),quantile(PhyCae.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.PhyCae.data)
List of 28
$ cov.dis : num [1:404] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:404] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.dis : num [1:404] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.dis : num [1:404] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.dis : num [1:404] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.dis : num [1:404] 0.796 1.722 1.028 1.375 0.796 ...
$ N_discrete : int 404
$ cov.cont : num [1:10] 0.04 0.08 0.08 0.12 0.08 0.04 0.04 0.12 0.44 0.12
$ plotgroup.cont: num [1:10] 9 11 14 17 24 30 36 51 51 69
$ sri.cont : num [1:10] -0.7393 -0.0933 -0.2978 0.4817 0.075 ...
$ tri.cont : num [1:10] -1.127 -0.199 0.634 -0.245 0.119 ...
$ twi.cont : num [1:10] -0.876 -0.637 2.047 -0.85 0.574 ...
$ compet.cont : num [1:10] -0.708 3.227 0.218 2.07 0.681 ...
$ N_cont : int 10
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.656 -0.603 -0.55 -0.498 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# RhoGro ----
shrub_gradient_jags.RhoGro.data <- list(
# plot level predictors, for discrete...
cov.dis = RhoGro.dis$cover,
plotgroup.dis = RhoGro.dis$plotgroup.NUM, #AB added this
# isocline.dis = RhoGro.dis$site_alt.NUM,
# inclin_down.dis = RhoGro.dis$inclin_downC,
sri.dis = RhoGro.dis$sriC,
tri.dis = RhoGro.dis$triC,
twi.dis = RhoGro.dis$twiC,
compet.dis = RhoGro.dis$competC,
N_discrete = nrow(RhoGro.dis),
# ...and continuous part of the data
cov.cont = RhoGro.cont$cover,
plotgroup.cont = RhoGro.cont$plotgroup.NUM, #AB added this
# isocline.cont = RhoGro.cont$site_alt.NUM,
# inclin_down.cont = RhoGro.cont$inclin_downC,
sri.cont = RhoGro.cont$sriC,
tri.cont = RhoGro.cont$triC,
twi.cont = RhoGro.cont$twiC,
compet.cont = RhoGro.cont$competC,
N_cont = nrow(RhoGro.cont),
# plot group level predictors
tempjja.tot = RhoGro.tot$tempjja_ts_30C[!duplicated(RhoGro.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = RhoGro.tot$tempmax_ts_30C[!duplicated(RhoGro.tot$plotgroup.NUM)],
# tempmin.tot = RhoGro.tot$tempmin_ts_30C[!duplicated(RhoGro.tot$plotgroup.NUM)],
tempcont.tot = RhoGro.tot$tempcont_ts_30C[!duplicated(RhoGro.tot$plotgroup.NUM)],
precipjja.tot = RhoGro.tot$precipjja_ts_30C[!duplicated(RhoGro.tot$plotgroup.NUM)],
# precipjfmam.tot = RhoGro.tot$precipjfmam_ts_30C[!duplicated(RhoGro.tot$plotgroup.NUM)]
N_plotgroups = length(unique(RhoGro.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = RhoGro.tot$altC[!duplicated(RhoGro.tot$site_alt.NUM)],
# N_isoclines = length(unique(RhoGro.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(RhoGro.tot$competC), to = max(RhoGro.tot$competC), length.out = 100),
xhat_sri = seq(from = min(RhoGro.tot$sriC), to = max(RhoGro.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(RhoGro.tot$triC), to = max(RhoGro.tot$triC), length.out = 100),
xhat_twi = seq(from = min(RhoGro.tot$twiC), to = max(RhoGro.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(RhoGro.tot$tempjja_ts_30C), to = max(RhoGro.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(RhoGro.tot$precipjja_ts_30C), to = max(RhoGro.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(RhoGro.tot$tempcont_ts_30C), to = max(RhoGro.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(RhoGro.tot$tempjja_ts_30C,0.05),quantile(RhoGro.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.RhoGro.data)
List of 28
$ cov.dis : num [1:325] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:325] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.dis : num [1:325] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.dis : num [1:325] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.dis : num [1:325] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.dis : num [1:325] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_discrete : int 325
$ cov.cont : num [1:89] 0.08 0.12 0.52 0.36 0.28 0.04 0.04 0.16 0.04 0.04 ...
$ plotgroup.cont: num [1:89] 5 11 20 20 20 20 20 20 24 24 ...
$ sri.cont : num [1:89] 0.8231 -0.9499 -0.0661 -0.0445 -2.5788 ...
$ tri.cont : num [1:89] 0.4613 -0.0943 1.3399 1.4254 0.9241 ...
$ twi.cont : num [1:89] 1.775 -0.637 -1.027 -1.119 -1.119 ...
$ compet.cont : num [1:89] -0.708 -0.361 -0.708 -0.708 -0.708 ...
$ N_cont : int 89
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.697 -0.685 -0.673 -0.661 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# RhoTom ----
shrub_gradient_jags.RhoTom.data <- list(
# plot level predictors, for discrete...
cov.dis = RhoTom.dis$cover,
plotgroup.dis = RhoTom.dis$plotgroup.NUM, #AB added this
# isocline.dis = RhoTom.dis$site_alt.NUM,
# inclin_down.dis = RhoTom.dis$inclin_downC,
sri.dis = RhoTom.dis$sriC,
tri.dis = RhoTom.dis$triC,
twi.dis = RhoTom.dis$twiC,
compet.dis = RhoTom.dis$competC,
N_discrete = nrow(RhoTom.dis),
# ...and continuous part of the data
cov.cont = RhoTom.cont$cover,
plotgroup.cont = RhoTom.cont$plotgroup.NUM, #AB added this
# isocline.cont = RhoTom.cont$site_alt.NUM,
# inclin_down.cont = RhoTom.cont$inclin_downC,
sri.cont = RhoTom.cont$sriC,
tri.cont = RhoTom.cont$triC,
twi.cont = RhoTom.cont$twiC,
compet.cont = RhoTom.cont$competC,
N_cont = nrow(RhoTom.cont),
# plot group level predictors
tempjja.tot = RhoTom.tot$tempjja_ts_30C[!duplicated(RhoTom.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = RhoTom.tot$tempmax_ts_30C[!duplicated(RhoTom.tot$plotgroup.NUM)],
# tempmin.tot = RhoTom.tot$tempmin_ts_30C[!duplicated(RhoTom.tot$plotgroup.NUM)],
tempcont.tot = RhoTom.tot$tempcont_ts_30C[!duplicated(RhoTom.tot$plotgroup.NUM)],
precipjja.tot = RhoTom.tot$precipjja_ts_30C[!duplicated(RhoTom.tot$plotgroup.NUM)],
# precipjfmam.tot = RhoTom.tot$precipjfmam_ts_30C[!duplicated(RhoTom.tot$plotgroup.NUM)]
N_plotgroups = length(unique(RhoTom.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = RhoTom.tot$altC[!duplicated(RhoTom.tot$site_alt.NUM)],
# N_isoclines = length(unique(RhoTom.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(RhoTom.tot$competC), to = max(RhoTom.tot$competC), length.out = 100),
xhat_sri = seq(from = min(RhoTom.tot$sriC), to = max(RhoTom.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(RhoTom.tot$triC), to = max(RhoTom.tot$triC), length.out = 100),
xhat_twi = seq(from = min(RhoTom.tot$twiC), to = max(RhoTom.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(RhoTom.tot$tempjja_ts_30C), to = max(RhoTom.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(RhoTom.tot$precipjja_ts_30C), to = max(RhoTom.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(RhoTom.tot$tempcont_ts_30C), to = max(RhoTom.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(RhoTom.tot$tempjja_ts_30C,0.05),quantile(RhoTom.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.RhoTom.data)
List of 28
$ cov.dis : num [1:401] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:401] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.dis : num [1:401] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.dis : num [1:401] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.dis : num [1:401] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.dis : num [1:401] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_discrete : int 401
$ cov.cont : num [1:13] 0.12 0.2 0.12 0.4 0.2 0.2 0.2 0.24 0.04 0.2 ...
$ plotgroup.cont: num [1:13] 11 11 12 30 30 36 39 39 42 43 ...
$ sri.cont : num [1:13] -0.95 0.723 -0.185 0.633 -1.308 ...
$ tri.cont : num [1:13] -0.0943 -0.2387 -0.4298 1.7289 1.5645 ...
$ twi.cont : num [1:13] -0.637 -0.637 0.549 -1.385 -1.362 ...
$ compet.cont : num [1:13] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_cont : int 13
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# SalArc ----
shrub_gradient_jags.SalArc.data <- list(
# plot level predictors, for discrete...
cov.dis = SalArc.dis$cover,
plotgroup.dis = SalArc.dis$plotgroup.NUM, #AB added this
# isocline.dis = SalArc.dis$site_alt.NUM,
# inclin_down.dis = SalArc.dis$inclin_downC,
sri.dis = SalArc.dis$sriC,
tri.dis = SalArc.dis$triC,
twi.dis = SalArc.dis$twiC,
compet.dis = SalArc.dis$competC,
N_discrete = nrow(SalArc.dis),
# ...and continuous part of the data
cov.cont = SalArc.cont$cover,
plotgroup.cont = SalArc.cont$plotgroup.NUM, #AB added this
# isocline.cont = SalArc.cont$site_alt.NUM,
# inclin_down.cont = SalArc.cont$inclin_downC,
sri.cont = SalArc.cont$sriC,
tri.cont = SalArc.cont$triC,
twi.cont = SalArc.cont$twiC,
compet.cont = SalArc.cont$competC,
N_cont = nrow(SalArc.cont),
# plot group level predictors
tempjja.tot = SalArc.tot$tempjja_ts_30C[!duplicated(SalArc.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = SalArc.tot$tempmax_ts_30C[!duplicated(SalArc.tot$plotgroup.NUM)],
# tempmin.tot = SalArc.tot$tempmin_ts_30C[!duplicated(SalArc.tot$plotgroup.NUM)],
tempcont.tot = SalArc.tot$tempcont_ts_30C[!duplicated(SalArc.tot$plotgroup.NUM)],
precipjja.tot = SalArc.tot$precipjja_ts_30C[!duplicated(SalArc.tot$plotgroup.NUM)],
# precipjfmam.tot = SalArc.tot$precipjfmam_ts_30C[!duplicated(SalArc.tot$plotgroup.NUM)]
N_plotgroups = length(unique(SalArc.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = SalArc.tot$altC[!duplicated(SalArc.tot$site_alt.NUM)],
# N_isoclines = length(unique(SalArc.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(SalArc.tot$competC), to = max(SalArc.tot$competC), length.out = 100),
xhat_sri = seq(from = min(SalArc.tot$sriC), to = max(SalArc.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(SalArc.tot$triC), to = max(SalArc.tot$triC), length.out = 100),
xhat_twi = seq(from = min(SalArc.tot$twiC), to = max(SalArc.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(SalArc.tot$tempjja_ts_30C), to = max(SalArc.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(SalArc.tot$precipjja_ts_30C), to = max(SalArc.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(SalArc.tot$tempcont_ts_30C), to = max(SalArc.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(SalArc.tot$tempjja_ts_30C,0.05),quantile(SalArc.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.SalArc.data)
List of 28
$ cov.dis : num [1:401] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:401] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.dis : num [1:401] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.dis : num [1:401] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.dis : num [1:401] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.dis : num [1:401] 0.796 1.722 1.028 1.375 0.796 ...
$ N_discrete : int 401
$ cov.cont : num [1:13] 0.04 0.04 0.08 0.04 0.04 0.04 0.04 0.2 0.08 0.04 ...
$ plotgroup.cont: num [1:13] 4 4 5 9 16 17 18 30 43 45 ...
$ sri.cont : num [1:13] -1.344 -2.825 0.823 0.14 0.305 ...
$ tri.cont : num [1:13] 1.498 1.327 0.461 -0.478 -0.95 ...
$ twi.cont : num [1:13] 1.132 1.132 1.775 -0.876 -1.179 ...
$ compet.cont : num [1:13] 1.028 0.796 2.88 0.565 -0.13 ...
$ N_cont : int 13
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.656 -0.603 -0.55 -0.498 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# SalGla ----
shrub_gradient_jags.SalGla.data <- list(
# plot level predictors, for discrete...
cov.dis = SalGla.dis$cover,
plotgroup.dis = SalGla.dis$plotgroup.NUM, #AB added this
# isocline.dis = SalGla.dis$site_alt.NUM,
# inclin_down.dis = SalGla.dis$inclin_downC,
sri.dis = SalGla.dis$sriC,
tri.dis = SalGla.dis$triC,
twi.dis = SalGla.dis$twiC,
compet.dis = SalGla.dis$competC,
N_discrete = nrow(SalGla.dis),
# ...and continuous part of the data
cov.cont = SalGla.cont$cover,
plotgroup.cont = SalGla.cont$plotgroup.NUM, #AB added this
# isocline.cont = SalGla.cont$site_alt.NUM,
# inclin_down.cont = SalGla.cont$inclin_downC,
sri.cont = SalGla.cont$sriC,
tri.cont = SalGla.cont$triC,
twi.cont = SalGla.cont$twiC,
compet.cont = SalGla.cont$competC,
N_cont = nrow(SalGla.cont),
# plot group level predictors
tempjja.tot = SalGla.tot$tempjja_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = SalGla.tot$tempmax_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
# tempmin.tot = SalGla.tot$tempmin_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
tempcont.tot = SalGla.tot$tempcont_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
precipjja.tot = SalGla.tot$precipjja_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
# precipjfmam.tot = SalGla.tot$precipjfmam_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)]
N_plotgroups = length(unique(SalGla.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = SalGla.tot$altC[!duplicated(SalGla.tot$site_alt.NUM)],
# N_isoclines = length(unique(SalGla.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(SalGla.tot$competC), to = max(SalGla.tot$competC), length.out = 100),
xhat_sri = seq(from = min(SalGla.tot$sriC), to = max(SalGla.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(SalGla.tot$triC), to = max(SalGla.tot$triC), length.out = 100),
xhat_twi = seq(from = min(SalGla.tot$twiC), to = max(SalGla.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(SalGla.tot$tempjja_ts_30C), to = max(SalGla.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(SalGla.tot$precipjja_ts_30C), to = max(SalGla.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(SalGla.tot$tempcont_ts_30C), to = max(SalGla.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(SalGla.tot$tempjja_ts_30C,0.05),quantile(SalGla.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.SalGla.data)
List of 28
$ cov.dis : num [1:289] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:289] 1 1 1 1 1 1 2 2 2 2 ...
$ sri.dis : num [1:289] 0.51 0.829 1.034 0.37 0.573 ...
$ tri.dis : num [1:289] -0.274 -0.318 -0.301 -0.274 -0.176 ...
$ twi.dis : num [1:289] 0.144 0.144 0.144 -0.384 -0.384 ...
$ compet.dis : num [1:289] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_discrete : int 289
$ cov.cont : num [1:125] 0.32 0.12 0.2 0.08 0.44 0.08 0.04 0.12 0.04 0.04 ...
$ plotgroup.cont: num [1:125] 2 2 4 4 6 8 11 11 12 14 ...
$ sri.cont : num [1:125] -1.562 -0.758 0.646 -2.68 -0.409 ...
$ tri.cont : num [1:125] -0.0114 0.1046 -0.9626 1.4106 0.1025 ...
$ twi.cont : num [1:125] -0.73 -0.73 1.13 1.13 1.96 ...
$ compet.cont : num [1:125] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_cont : int 125
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.686 -0.664 -0.642 -0.619 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
# VacUli ----
shrub_gradient_jags.VacUli.data <- list(
# plot level predictors, for discrete...
cov.dis = VacUli.dis$cover,
plotgroup.dis = VacUli.dis$plotgroup.NUM, #AB added this
# isocline.dis = VacUli.dis$site_alt.NUM,
# inclin_down.dis = VacUli.dis$inclin_downC,
sri.dis = VacUli.dis$sriC,
tri.dis = VacUli.dis$triC,
twi.dis = VacUli.dis$twiC,
compet.dis = VacUli.dis$competC,
N_discrete = nrow(VacUli.dis),
# ...and continuous part of the data
cov.cont = VacUli.cont$cover,
plotgroup.cont = VacUli.cont$plotgroup.NUM, #AB added this
# isocline.cont = VacUli.cont$site_alt.NUM,
# inclin_down.cont = VacUli.cont$inclin_downC,
sri.cont = VacUli.cont$sriC,
tri.cont = VacUli.cont$triC,
twi.cont = VacUli.cont$twiC,
compet.cont = VacUli.cont$competC,
N_cont = nrow(VacUli.cont),
# plot group level predictors
tempjja.tot = VacUli.tot$tempjja_ts_30C[!duplicated(VacUli.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = VacUli.tot$tempmax_ts_30C[!duplicated(VacUli.tot$plotgroup.NUM)],
# tempmin.tot = VacUli.tot$tempmin_ts_30C[!duplicated(VacUli.tot$plotgroup.NUM)],
tempcont.tot = VacUli.tot$tempcont_ts_30C[!duplicated(VacUli.tot$plotgroup.NUM)],
precipjja.tot = VacUli.tot$precipjja_ts_30C[!duplicated(VacUli.tot$plotgroup.NUM)],
# precipjfmam.tot = VacUli.tot$precipjfmam_ts_30C[!duplicated(VacUli.tot$plotgroup.NUM)]
N_plotgroups = length(unique(VacUli.tot$site_alt_plotgroup_id)),
# # site/alt level predictors
# alt.tot = VacUli.tot$altC[!duplicated(VacUli.tot$site_alt.NUM)],
# N_isoclines = length(unique(VacUli.tot$site_alt_id))
# subset of values for prediction, for each predictor...
xhat_compet = seq(from = min(VacUli.tot$competC), to = max(VacUli.tot$competC), length.out = 100),
xhat_sri = seq(from = min(VacUli.tot$sriC), to = max(VacUli.tot$sriC), length.out = 100),
xhat_tri = seq(from = min(VacUli.tot$triC), to = max(VacUli.tot$triC), length.out = 100),
xhat_twi = seq(from = min(VacUli.tot$twiC), to = max(VacUli.tot$twiC), length.out = 100),
xhat_tempjja = seq(from = min(VacUli.tot$tempjja_ts_30C), to = max(VacUli.tot$tempjja_ts_30C), length.out = 100),
xhat_precipjja = seq(from = min(VacUli.tot$precipjja_ts_30C), to = max(VacUli.tot$precipjja_ts_30C), length.out = 100),
xhat_tempcont = seq(from = min(VacUli.tot$tempcont_ts_30C), to = max(VacUli.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100,
# ... and for predicting at high/low temperature levels
xhat_tempjja2 = as.numeric(c(quantile(VacUli.tot$tempjja_ts_30C,0.05),quantile(VacUli.tot$tempjja_ts_30C,0.95))),
Nxhat2 = 2
)
str(shrub_gradient_jags.VacUli.data)
List of 28
$ cov.dis : num [1:309] 0 0 0 0 0 0 0 0 0 0 ...
$ plotgroup.dis : num [1:309] 1 1 1 1 2 2 2 2 2 2 ...
$ sri.dis : num [1:309] 0.51 0.829 1.034 0.573 0.945 ...
$ tri.dis : num [1:309] -0.274 -0.318 -0.301 -0.176 -0.983 ...
$ twi.dis : num [1:309] 0.144 0.144 0.144 -0.384 -0.754 ...
$ compet.dis : num [1:309] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_discrete : int 309
$ cov.cont : num [1:105] 0.08 0.04 0.24 0.28 0.08 0.32 0.36 0.28 0.08 0.16 ...
$ plotgroup.cont: num [1:105] 1 1 3 3 4 4 5 10 10 11 ...
$ sri.cont : num [1:105] 0.37 0.688 -0.207 -1.955 -0.556 ...
$ tri.cont : num [1:105] -0.274 0.187 0.138 0.355 1.316 ...
$ twi.cont : num [1:105] -0.384 -0.384 -1.069 -1.069 1.132 ...
$ compet.cont : num [1:105] -0.708 -0.708 -0.708 -0.708 -0.708 ...
$ N_cont : int 105
$ tempjja.tot : num [1:69] -0.321 -0.283 -0.226 -0.395 -0.42 ...
$ tempcont.tot : num [1:69] -2.45 -2.52 -2.49 -2.41 -2.39 ...
$ precipjja.tot : num [1:69] 0.713 0.866 1.044 0.626 0.593 ...
$ N_plotgroups : int 69
$ xhat_compet : num [1:100] -0.708 -0.667 -0.626 -0.585 -0.545 ...
$ xhat_sri : num [1:100] -3.1 -3.06 -3.02 -2.98 -2.94 ...
$ xhat_tri : num [1:100] -1.71 -1.65 -1.59 -1.53 -1.47 ...
$ xhat_twi : num [1:100] -1.7 -1.66 -1.62 -1.57 -1.53 ...
$ xhat_tempjja : num [1:100] -2.57 -2.52 -2.48 -2.44 -2.4 ...
$ xhat_precipjja: num [1:100] -1.66 -1.62 -1.58 -1.55 -1.51 ...
$ xhat_tempcont : num [1:100] -2.52 -2.49 -2.45 -2.42 -2.38 ...
$ Nxhat : num 100
$ xhat_tempjja2 : num [1:2] -2.01 1.31
$ Nxhat2 : num 2
write("
model{
# priors
intercept ~ dnorm(0, 0.0001)
b.compet ~ dnorm(0, 0.0001)
b.sri ~ dnorm(0, 0.0001)
b.tri ~ dnorm(0, 0.0001)
b.twi ~ dnorm(0, 0.0001)
sigma.plotgroup ~ dunif(0,100)
tau.plotgroup <- 1/(sigma.plotgroup * sigma.plotgroup)
b.tempjja.x ~ dnorm(0, 0.001)
b.tempjja.x2 ~ dnorm(0, 0.001)
b.tempcont.x ~ dnorm(0, 0.001)
b.tempcont.x2 ~ dnorm(0, 0.001)
b.precipjja.x ~ dnorm(0, 0.001)
b.precipjja.x2 ~ dnorm(0, 0.001)
phi ~ dgamma(0.1, 0.1)
# LIKELIHOOD for discrete part
for (i in 1:N_discrete){
cov.dis[i] ~ dbern(mu[i])
logit(mu[i]) <- b_plotgroup[plotgroup.dis[i]] + # ~= random effect of plot group
b.compet * compet.dis[i] +
b.twi * twi.dis[i] +
b.sri * sri.dis[i] +
b.tri * tri.dis[i]
}
# LIKELIHOOD for continuous part
for (j in 1:N_cont){
cov.cont[j] ~ dbeta(p[j], q[j])
p[j] <- mu2[j] * phi
q[j] <- (1 - mu2[j]) * phi
logit(mu2[j]) <- b_plotgroup[plotgroup.cont[j]] + # ~= random effect of plot group
b.compet * compet.cont[j] +
b.twi * twi.cont[j] +
b.sri * sri.cont[j] +
b.tri * tri.cont[j]
}
for (k in 1:N_plotgroups){ # length of total plotgroups
b_plotgroup[k] ~ dnorm(mu.plotgroup[k],tau.plotgroup)
mu.plotgroup[k] <- intercept +
# plot group level predictors, linear and quadratic term
b.tempjja.x * tempjja.tot[k] +
b.tempjja.x2 * (tempjja.tot[k]^2) +
b.tempcont.x * tempcont.tot[k] +
b.tempcont.x2 * (tempcont.tot[k]^2) +
b.precipjja.x * precipjja.tot[k] +
b.precipjja.x2 * (precipjja.tot[k]^2)
}
# add predicted values (derived parameters)
for (m in 1:Nxhat){
phat_compet[m] <- intercept + b.compet * xhat_compet[m]
phat_sri[m] <- intercept + b.sri * xhat_sri[m]
phat_tri[m] <- intercept + b.tri * xhat_tri[m]
phat_twi[m] <- intercept + b.twi * xhat_twi[m]
phat_tempjja[m] <- intercept + b.tempjja.x * xhat_tempjja[m] + b.tempjja.x2 * (xhat_tempjja[m]^2)
phat_tempcont[m] <- intercept + b.tempcont.x * xhat_tempcont[m] + b.tempcont.x2 * (xhat_tempcont[m]^2)
phat_precipjja[m] <- intercept + b.precipjja.x * xhat_precipjja[m] + b.precipjja.x2 * (xhat_precipjja[m]^2)
}
}
", file.path("..", "models", "shrub_gradient.spec.jags"))
Specify the parameters to be monitored:
params <- c("intercept",
"b.tempjja.x", "b.tempjja.x2",
"b.tempcont.x", "b.tempcont.x2",
"b.precipjja.x", "b.precipjja.x2",
"b.compet",
"b.sri",
"b.tri",
"b.twi",
"b_plotgroup[1]","b_plotgroup[2]","b_plotgroup[3]","b_plotgroup[63]",
"sigma.plotgroup",
"phi",
"phat_compet", "phat_sri", "phat_tri", "phat_twi", "phat_tempjja", "phat_tempcont", "phat_precipjja")
# run model
model_out.shrub_gradient.BetNan <- jags(shrub_gradient_jags.BetNan.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 8536
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.BetNan) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.BetNan <- model_out.shrub_gradient.BetNan$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.BetNan),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.BetNan <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.BetNan$BUGSoutput$sims.list)-4)){
ci_90.BetNan[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.BetNan$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.BetNan[param, 3] <- names(data.frame(model_out.shrub_gradient.BetNan$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.BetNan <- coeff.shrub_gradient.BetNan %>%
left_join(ci_90.BetNan, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
(effect_size_plot.BetNan <- model_plot_marg_function(coeff.shrub_gradient.BetNan, title_string = "Betula nana", plot_width = 11))
As the quadratic terms for summer temperature and precipitation were not significant, they were removed from the model before re-running:
# new model object with terms removed
write("
model{
# priors
intercept ~ dnorm(0, 0.0001)
b.compet ~ dnorm(0, 0.0001)
b.sri ~ dnorm(0, 0.0001)
b.tri ~ dnorm(0, 0.0001)
b.twi ~ dnorm(0, 0.0001)
sigma.plotgroup ~ dunif(0,100)
tau.plotgroup <- 1/(sigma.plotgroup * sigma.plotgroup)
b.tempjja.x ~ dnorm(0, 0.001)
# b.tempjja.x2 ~ dnorm(0, 0.001)
b.tempcont.x ~ dnorm(0, 0.001)
b.tempcont.x2 ~ dnorm(0, 0.001)
b.precipjja.x ~ dnorm(0, 0.001)
# b.precipjja.x2 ~ dnorm(0, 0.001)
phi ~ dgamma(0.1, 0.1)
# LIKELIHOOD for discrete part
for (i in 1:N_discrete){
cov.dis[i] ~ dbern(mu[i])
logit(mu[i]) <- b_plotgroup[plotgroup.dis[i]] + # ~= random effect of plot group
b.compet * compet.dis[i] +
b.twi * twi.dis[i] +
b.sri * sri.dis[i] +
b.tri * tri.dis[i]
}
# LIKELIHOOD for continuous part
for (j in 1:N_cont){
cov.cont[j] ~ dbeta(p[j], q[j])
p[j] <- mu2[j] * phi
q[j] <- (1 - mu2[j]) * phi
logit(mu2[j]) <- b_plotgroup[plotgroup.cont[j]] + # ~= random effect of plot group
b.compet * compet.cont[j] +
b.twi * twi.cont[j] +
b.sri * sri.cont[j] +
b.tri * tri.cont[j]
}
for (k in 1:N_plotgroups){ # length of total plotgroups
b_plotgroup[k] ~ dnorm(mu.plotgroup[k],tau.plotgroup)
mu.plotgroup[k] <- intercept +
# plot group level predictors, linear and quadratic term
b.tempjja.x * tempjja.tot[k] +
# b.tempjja.x2 * (tempjja.tot[k]^2) +
b.tempcont.x * tempcont.tot[k] +
b.tempcont.x2 * (tempcont.tot[k]^2) +
b.precipjja.x * precipjja.tot[k] # +
# b.precipjja.x2 * (precipjja.tot[k]^2)
}
# add predicted values (derived parameters)
for (m in 1:Nxhat){
phat_compet[m] <- intercept + b.compet * xhat_compet[m]
phat_sri[m] <- intercept +
b.sri * xhat_sri[m]
phat_tri[m] <- intercept +
b.tri * xhat_tri[m]
phat_twi[m] <- intercept +
b.twi * xhat_twi[m]
phat_tempjja[m] <- intercept +
b.tempjja.x * xhat_tempjja[m] # +
# b.tempjja.x2 * (xhat_tempjja[m]^2)
phat_tempcont[m] <- intercept +
b.tempcont.x * xhat_tempcont[m] +
b.tempcont.x2 * (xhat_tempcont[m]^2)
phat_precipjja[m] <- intercept +
b.precipjja.x * xhat_precipjja[m] # +
# b.precipjja.x2 * (xhat_precipjja[m]^2)
for (p in 1:Nxhat2){
phat_tempcontXtemp[m,p] <- intercept +
b.tempcont.x * xhat_tempcont[m] +
b.tempcont.x2 * (xhat_tempcont[m]^2) +
b.tempjja.x * xhat_tempjja2[p] # +
# b.tempjja.x2 * (xhat_tempjja2[p]^2)
}
}
}
", file.path("..", "models", "shrub_gradient.BetNan2.jags"))
# specify new set of parameters to be monitored
params_BetNan2 <- c("intercept",
"b.tempjja.x", # "b.tempjja.x2",
"b.tempcont.x", "b.tempcont.x2",
"b.precipjja.x", # "b.precipjja.x2",
"b.compet",
"b.sri",
"b.tri",
"b.twi",
"b_plotgroup[1]","b_plotgroup[2]","b_plotgroup[3]","b_plotgroup[63]",
"sigma.plotgroup",
"phi",
"phat_compet", "phat_sri", "phat_tri", "phat_twi", "phat_tempjja", "phat_tempcont", "phat_precipjja", "phat_tempcontXtemp")
model_out.shrub_gradient.BetNan2 <- jags(shrub_gradient_jags.BetNan.data, # input data
inits = NULL, # JAGS to create initial values
params_BetNan2, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.BetNan2.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 80
Total graph size: 8065
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.BetNan2) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.BetNan2 <- model_out.shrub_gradient.BetNan2$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# add 90% CIs
ci_90.BetNan2 <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.BetNan2$BUGSoutput$sims.list)-4)){
ci_90.BetNan2[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.BetNan2$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.BetNan2[param, 3] <- names(data.frame(model_out.shrub_gradient.BetNan2$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.BetNan2 <- coeff.shrub_gradient.BetNan2 %>%
left_join(ci_90.BetNan2, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.BetNan2 <- model_plot_marg_function(coeff.shrub_gradient.BetNan2, title_string = "Betula nana", plot_width = 8.5))
model_out.shrub_gradient.CasTet <- jags(shrub_gradient_jags.CasTet.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 8016
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.CasTet) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.CasTet <- model_out.shrub_gradient.CasTet$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.CasTet),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.CasTet <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.CasTet$BUGSoutput$sims.list)-4)){
ci_90.CasTet[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.CasTet$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.CasTet[param, 3] <- names(data.frame(model_out.shrub_gradient.CasTet$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.CasTet <- coeff.shrub_gradient.CasTet %>%
left_join(ci_90.CasTet, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# (effect_size_plot.CasTet <- model_plot_function(coeff.shrub_gradient.CasTet))
model_out.shrub_gradient.EmpNig <- jags(shrub_gradient_jags.EmpNig.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 8528
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.EmpNig) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.EmpNig <- model_out.shrub_gradient.EmpNig$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.EmpNig),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.EmpNig <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.EmpNig$BUGSoutput$sims.list)-4)){
ci_90.EmpNig[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.EmpNig$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.EmpNig[param, 3] <- names(data.frame(model_out.shrub_gradient.EmpNig$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.EmpNig <- coeff.shrub_gradient.EmpNig %>%
left_join(ci_90.EmpNig, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.EmpNig <- model_plot_sig_function(coeff.shrub_gradient.EmpNig, title_string = "Empetrum nigrum", plot_width = 11))
As the quadratic terms for temperature variability and precipitation were not significant, they were removed from the model before re-running:
# new model object with terms removed
write("
model{
# priors
intercept ~ dnorm(0, 0.0001)
b.compet ~ dnorm(0, 0.0001)
b.sri ~ dnorm(0, 0.0001)
b.tri ~ dnorm(0, 0.0001)
b.twi ~ dnorm(0, 0.0001)
sigma.plotgroup ~ dunif(0,100)
tau.plotgroup <- 1/(sigma.plotgroup * sigma.plotgroup)
b.tempjja.x ~ dnorm(0, 0.001)
b.tempjja.x2 ~ dnorm(0, 0.001)
b.tempcont.x ~ dnorm(0, 0.001)
# b.tempcont.x2 ~ dnorm(0, 0.001)
b.precipjja.x ~ dnorm(0, 0.001)
# b.precipjja.x2 ~ dnorm(0, 0.001)
phi ~ dgamma(0.1, 0.1)
# LIKELIHOOD for discrete part
for (i in 1:N_discrete){
cov.dis[i] ~ dbern(mu[i])
logit(mu[i]) <- b_plotgroup[plotgroup.dis[i]] + #AB added this, ~= random effect of plot group
b.compet * compet.dis[i] +
b.twi * twi.dis[i] +
b.sri * sri.dis[i] +
b.tri * tri.dis[i]
}
# LIKELIHOOD for continuous part
for (j in 1:N_cont){
cov.cont[j] ~ dbeta(p[j], q[j])
p[j] <- mu2[j] * phi
q[j] <- (1 - mu2[j]) * phi
logit(mu2[j]) <- b_plotgroup[plotgroup.cont[j]] + #AB added this, ~= random effect of plot group
b.compet * compet.cont[j] +
b.twi * twi.cont[j] +
b.sri * sri.cont[j] +
b.tri * tri.cont[j]
}
for (k in 1:N_plotgroups){ # length of total plotgroups
b_plotgroup[k] ~ dnorm(mu.plotgroup[k],tau.plotgroup)
mu.plotgroup[k] <- intercept +
# plot group level predictors, linear and quadratic term
b.tempjja.x * tempjja.tot[k] +
b.tempjja.x2 * (tempjja.tot[k]^2) +
b.tempcont.x * tempcont.tot[k] +
# b.tempcont.x2 * (tempcont.tot[k]^2) +
b.precipjja.x * precipjja.tot[k] # +
# b.precipjja.x2 * (precipjja.tot[k]^2)
}
# add predicted values (derived parameters)
for (m in 1:Nxhat){
phat_compet[m] <- intercept + b.compet * xhat_compet[m]
phat_sri[m] <- intercept +
b.sri * xhat_sri[m]
phat_tri[m] <- intercept +
b.tri * xhat_tri[m]
phat_twi[m] <- intercept +
b.twi * xhat_twi[m]
phat_tempjja[m] <- intercept +
b.tempjja.x * xhat_tempjja[m] +
b.tempjja.x2 * (xhat_tempjja[m]^2)
phat_tempcont[m] <- intercept +
b.tempcont.x * xhat_tempcont[m] # +
# b.tempcont.x2 * (xhat_tempcont[m]^2)
phat_precipjja[m] <- intercept +
b.precipjja.x * xhat_precipjja[m] # +
# b.precipjja.x2 * (xhat_precipjja[m]^2)
for (p in 1:Nxhat2){
phat_tempcontXtemp[m,p] <- intercept +
b.tempcont.x * xhat_tempcont[m] +
# b.tempcont.x2 * (xhat_tempcont[m]^2) +
b.tempjja.x * xhat_tempjja2[p] +
b.tempjja.x2 * (xhat_tempjja2[p]^2)
}
}
}
", file.path("..", "models", "shrub_gradient.EmpNig2.jags"))
# specify new set of parameters to be monitored
params_EmpNig2 <- c("intercept",
"b.tempjja.x", "b.tempjja.x2",
"b.tempcont.x", # "b.tempcont.x2",
"b.precipjja.x", # "b.precipjja.x2",
"b.compet",
"b.sri",
"b.tri",
"b.twi",
"b_plotgroup[1]","b_plotgroup[2]","b_plotgroup[3]","b_plotgroup[63]",
"sigma.plotgroup",
"phi",
"phat_compet", "phat_sri", "phat_tri", "phat_twi", "phat_tempjja", "phat_tempcont", "phat_precipjja", "phat_tempcontXtemp")
model_out.shrub_gradient.EmpNig2 <- jags(shrub_gradient_jags.EmpNig.data, # input data
inits = NULL, # JAGS to create initial values
params_EmpNig2, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.EmpNig2.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 80
Total graph size: 8063
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.EmpNig2) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.EmpNig2 <- model_out.shrub_gradient.EmpNig2$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# add 90% CIs
ci_90.EmpNig2 <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.EmpNig$BUGSoutput$sims.list)-4)){
ci_90.EmpNig2[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.EmpNig2$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.EmpNig2[param, 3] <- names(data.frame(model_out.shrub_gradient.EmpNig2$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.EmpNig2 <- coeff.shrub_gradient.EmpNig2 %>%
left_join(ci_90.EmpNig2, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.EmpNig2 <- model_plot_marg_function(coeff.shrub_gradient.EmpNig2, title_string = "Empetrum nigrum", plot_width = 8.5))
model_out.shrub_gradient.PhyCae <- jags(shrub_gradient_jags.PhyCae.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 7999
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.PhyCae) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.PhyCae <- model_out.shrub_gradient.PhyCae$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.PhyCae),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.PhyCae <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.PhyCae$BUGSoutput$sims.list)-4)){
ci_90.PhyCae[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.PhyCae$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.PhyCae[param, 3] <- names(data.frame(model_out.shrub_gradient.PhyCae$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.PhyCae <- coeff.shrub_gradient.PhyCae %>%
left_join(ci_90.PhyCae, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
NA
model_out.shrub_gradient.RhoGro <- jags(shrub_gradient_jags.RhoGro.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 8210
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.RhoGro) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.RhoGro <- model_out.shrub_gradient.RhoGro$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.RhoGro),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.RhoGro <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.RhoGro$BUGSoutput$sims.list)-4)){
ci_90.RhoGro[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.RhoGro$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.RhoGro[param, 3] <- names(data.frame(model_out.shrub_gradient.RhoGro$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.RhoGro <- coeff.shrub_gradient.RhoGro %>%
left_join(ci_90.RhoGro, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.RhoGro <- model_plot_sig_function(coeff.shrub_gradient.RhoGro, title_string = "Rhododendron groenlandicum", plot_width = 11))
As the quadratic terms for temperature variability and precipitation were not significant, they were removed from the model before re-running:
# new model object with terms removed
write("
model{
# priors
intercept ~ dnorm(0, 0.0001)
b.compet ~ dnorm(0, 0.0001)
b.sri ~ dnorm(0, 0.0001)
b.tri ~ dnorm(0, 0.0001)
b.twi ~ dnorm(0, 0.0001)
sigma.plotgroup ~ dunif(0,100)
tau.plotgroup <- 1/(sigma.plotgroup * sigma.plotgroup)
b.tempjja.x ~ dnorm(0, 0.001)
b.tempjja.x2 ~ dnorm(0, 0.001)
b.tempcont.x ~ dnorm(0, 0.001)
# b.tempcont.x2 ~ dnorm(0, 0.001)
b.precipjja.x ~ dnorm(0, 0.001)
# b.precipjja.x2 ~ dnorm(0, 0.001)
phi ~ dgamma(0.1, 0.1)
# LIKELIHOOD for discrete part
for (i in 1:N_discrete){
cov.dis[i] ~ dbern(mu[i])
logit(mu[i]) <- b_plotgroup[plotgroup.dis[i]] + # ~= random effect of plot group
b.compet * compet.dis[i] +
b.twi * twi.dis[i] +
b.sri * sri.dis[i] +
b.tri * tri.dis[i]
}
# LIKELIHOOD for continuous part
for (j in 1:N_cont){
cov.cont[j] ~ dbeta(p[j], q[j])
p[j] <- mu2[j] * phi
q[j] <- (1 - mu2[j]) * phi
logit(mu2[j]) <- b_plotgroup[plotgroup.cont[j]] + # ~= random effect of plot group
b.compet * compet.cont[j] +
b.twi * twi.cont[j] +
b.sri * sri.cont[j] +
b.tri * tri.cont[j]
}
for (k in 1:N_plotgroups){ # length of total plotgroups
b_plotgroup[k] ~ dnorm(mu.plotgroup[k],tau.plotgroup)
mu.plotgroup[k] <- intercept +
# plot group level predictors, linear and quadratic term
b.tempjja.x * tempjja.tot[k] +
b.tempjja.x2 * (tempjja.tot[k]^2) +
b.tempcont.x * tempcont.tot[k] +
# b.tempcont.x2 * (tempcont.tot[k]^2) +
b.precipjja.x * precipjja.tot[k] # +
# b.precipjja.x2 * (precipjja.tot[k]^2)
}
# add predicted values (derived parameters)
for (m in 1:Nxhat){
phat_compet[m] <- intercept + b.compet * xhat_compet[m]
phat_sri[m] <- intercept +
b.sri * xhat_sri[m]
phat_tri[m] <- intercept +
b.tri * xhat_tri[m]
phat_twi[m] <- intercept +
b.twi * xhat_twi[m]
phat_tempjja[m] <- intercept +
b.tempjja.x * xhat_tempjja[m] +
b.tempjja.x2 * (xhat_tempjja[m]^2)
phat_tempcont[m] <- intercept +
b.tempcont.x * xhat_tempcont[m] # +
# b.tempcont.x2 * (xhat_tempcont[m]^2)
phat_precipjja[m] <- intercept +
b.precipjja.x * xhat_precipjja[m] # +
# b.precipjja.x2 * (xhat_precipjja[m]^2)
for (p in 1:Nxhat2){
phat_tempcontXtemp[m,p] <- intercept +
b.tempcont.x * xhat_tempcont[m] +
# b.tempcont.x2 * (xhat_tempcont[m]^2) +
b.tempjja.x * xhat_tempjja2[p] +
b.tempjja.x2 * (xhat_tempjja2[p]^2)
}
}
}
", file.path("..", "models", "shrub_gradient.RhoGro2.jags"))
# specify new set of parameters to be monitored
params_RhoGro2 <- c("intercept",
"b.tempjja.x", "b.tempjja.x2",
"b.tempcont.x", # "b.tempcont.x2",
"b.precipjja.x", # "b.precipjja.x2",
"b.compet",
"b.sri",
"b.tri",
"b.twi",
"b_plotgroup[1]","b_plotgroup[2]","b_plotgroup[3]","b_plotgroup[63]",
"sigma.plotgroup",
"phi",
"phat_compet", "phat_sri", "phat_tri", "phat_twi", "phat_tempjja", "phat_tempcont", "phat_precipjja", "phat_tempcontXtemp")
model_out.shrub_gradient.RhoGro2 <- jags(shrub_gradient_jags.RhoGro.data, # input data
inits = NULL, # JAGS to create initial values
params_RhoGro2, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.RhoGro2.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 80
Total graph size: 7745
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.RhoGro2) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.RhoGro2 <- model_out.shrub_gradient.RhoGro2$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.RhoGro),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.RhoGro2 <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.RhoGro$BUGSoutput$sims.list)-4)){
ci_90.RhoGro2[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.RhoGro2$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.RhoGro2[param, 3] <- names(data.frame(model_out.shrub_gradient.RhoGro2$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.RhoGro2 <- coeff.shrub_gradient.RhoGro2 %>%
left_join(ci_90.RhoGro2, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.RhoGro2 <- model_plot_sig_function(coeff.shrub_gradient.RhoGro2, title_string = "Rhododendron groenlandicum", plot_width = 8.5))
model_out.shrub_gradient.RhoTom <- jags(shrub_gradient_jags.RhoTom.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 7780
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.RhoTom) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.RhoTom <- model_out.shrub_gradient.RhoTom$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.RhoTom),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.RhoTom <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.RhoTom$BUGSoutput$sims.list)-4)){
ci_90.RhoTom[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.RhoTom$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.RhoTom[param, 3] <- names(data.frame(model_out.shrub_gradient.RhoTom$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.RhoTom <- coeff.shrub_gradient.RhoTom %>%
left_join(ci_90.RhoTom, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
NA
model_out.shrub_gradient.SalArc <- jags(shrub_gradient_jags.SalArc.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 8008
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.SalArc) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.SalArc <- model_out.shrub_gradient.SalArc$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.SalArc),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.SalArc <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.SalArc$BUGSoutput$sims.list)-4)){
ci_90.SalArc[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.SalArc$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.SalArc[param, 3] <- names(data.frame(model_out.shrub_gradient.SalArc$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.SalArc <- coeff.shrub_gradient.SalArc %>%
left_join(ci_90.SalArc, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
NA
model_out.shrub_gradient.SalGla <- jags(shrub_gradient_jags.SalGla.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 8331
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.SalGla) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.SalGla <- model_out.shrub_gradient.SalGla$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.SalGla),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.SalGla <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.SalGla$BUGSoutput$sims.list)-4)){
ci_90.SalGla[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.SalGla$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.SalGla[param, 3] <- names(data.frame(model_out.shrub_gradient.SalGla$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.SalGla <- coeff.shrub_gradient.SalGla %>%
left_join(ci_90.SalGla, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.SalGla <- model_plot_sig_function(coeff.shrub_gradient.SalGla, title_string = "Salix glauca", plot_width = 11))
As the quadratic terms for summer temperature, temperature variability and summer precipitation were not significant, they were removed from the model before re-running:
# new model object with terms removed
write("
model{
# priors
intercept ~ dnorm(0, 0.0001)
b.compet ~ dnorm(0, 0.0001)
b.sri ~ dnorm(0, 0.0001)
b.tri ~ dnorm(0, 0.0001)
b.twi ~ dnorm(0, 0.0001)
sigma.plotgroup ~ dunif(0,100)
tau.plotgroup <- 1/(sigma.plotgroup * sigma.plotgroup)
b.tempjja.x ~ dnorm(0, 0.001)
b.tempjja.x2 ~ dnorm(0, 0.001)
b.tempcont.x ~ dnorm(0, 0.001)
# b.tempcont.x2 ~ dnorm(0, 0.001)
b.precipjja.x ~ dnorm(0, 0.001)
# b.precipjja.x2 ~ dnorm(0, 0.001)
phi ~ dgamma(0.1, 0.1)
# LIKELIHOOD for discrete part
for (i in 1:N_discrete){
cov.dis[i] ~ dbern(mu[i])
logit(mu[i]) <- b_plotgroup[plotgroup.dis[i]] + # ~= random effect of plot group
b.compet * compet.dis[i] +
b.twi * twi.dis[i] +
b.sri * sri.dis[i] +
b.tri * tri.dis[i]
}
# LIKELIHOOD for continuous part
for (j in 1:N_cont){
cov.cont[j] ~ dbeta(p[j], q[j])
p[j] <- mu2[j] * phi
q[j] <- (1 - mu2[j]) * phi
logit(mu2[j]) <- b_plotgroup[plotgroup.cont[j]] + # ~= random effect of plot group
b.compet * compet.cont[j] +
b.twi * twi.cont[j] +
b.sri * sri.cont[j] +
b.tri * tri.cont[j]
}
for (k in 1:N_plotgroups){ # length of total plotgroups
b_plotgroup[k] ~ dnorm(mu.plotgroup[k],tau.plotgroup)
mu.plotgroup[k] <- intercept +
# plot group level predictors, linear and quadratic term
b.tempjja.x * tempjja.tot[k] +
# b.tempjja.x2 * (tempjja.tot[k]^2) +
b.tempcont.x * tempcont.tot[k] +
# b.tempcont.x2 * (tempcont.tot[k]^2) +
b.precipjja.x * precipjja.tot[k] # +
# b.precipjja.x2 * (precipjja.tot[k]^2)
}
# add predicted values (derived parameters)
for (m in 1:Nxhat){
phat_compet[m] <- intercept + b.compet * xhat_compet[m]
phat_sri[m] <- intercept +
b.sri * xhat_sri[m]
phat_tri[m] <- intercept +
b.tri * xhat_tri[m]
phat_twi[m] <- intercept +
b.twi * xhat_twi[m]
phat_tempjja[m] <- intercept +
b.tempjja.x * xhat_tempjja[m] # +
# b.tempjja.x2 * (xhat_tempjja[m]^2)
phat_tempcont[m] <- intercept +
b.tempcont.x * xhat_tempcont[m] # +
# b.tempcont.x2 * (xhat_tempcont[m]^2)
phat_precipjja[m] <- intercept +
b.precipjja.x * xhat_precipjja[m] # +
# b.precipjja.x2 * (xhat_precipjja[m]^2)
for (p in 1:Nxhat2){
phat_tempcontXtemp[m,p] <- intercept +
b.tempcont.x * xhat_tempcont[m] +
# b.tempcont.x2 * (xhat_tempcont[m]^2) +
b.tempjja.x * xhat_tempjja2[p] # +
# b.tempjja.x2 * (xhat_tempjja2[p]^2)
}
}
}
", file.path("..", "models", "shrub_gradient.SalGla2.jags"))
# specify new set of parameters to be monitored
params_SalGla2 <- c("intercept",
"b.tempjja.x", # "b.tempjja.x2",
"b.tempcont.x", # "b.tempcont.x2",
"b.precipjja.x", # "b.precipjja.x2",
"b.compet",
"b.sri",
"b.tri",
"b.twi",
"b_plotgroup[1]","b_plotgroup[2]","b_plotgroup[3]","b_plotgroup[63]",
"sigma.plotgroup",
"phi",
"phat_compet", "phat_sri", "phat_tri", "phat_twi", "phat_tempjja", "phat_tempcont", "phat_precipjja", "phat_tempcontXtemp")
model_out.shrub_gradient.SalGla2 <- jags(shrub_gradient_jags.SalGla.data, # input data
inits = NULL, # JAGS to create initial values
params_SalGla2, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.SalGla2.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 80
Total graph size: 7523
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.SalGla2) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.SalGla2 <- model_out.shrub_gradient.SalGla2$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# add 90% CIs
ci_90.SalGla2 <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.SalGla$BUGSoutput$sims.list)-4)){
ci_90.SalGla2[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.SalGla2[param, 3] <- names(data.frame(model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.SalGla2 <- coeff.shrub_gradient.SalGla2 %>%
left_join(ci_90.SalGla2, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.SalGla2 <- model_plot_marg_function(coeff.shrub_gradient.SalGla2, title_string = "Salix glauca", plot_width = 8))
model_out.shrub_gradient.VacUli <- jags(shrub_gradient_jags.VacUli.data, # input data
inits = NULL, # JAGS to create initial values
params, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.spec.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Unused variable "xhat_tempjja2" in dataUnused variable "Nxhat2" in data
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 82
Total graph size: 8284
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.VacUli) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.VacUli <- model_out.shrub_gradient.VacUli$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.VacUli),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.VacUli <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.VacUli$BUGSoutput$sims.list)-4)){
ci_90.VacUli[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.VacUli$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.VacUli[param, 3] <- names(data.frame(model_out.shrub_gradient.VacUli$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.VacUli <- coeff.shrub_gradient.VacUli %>%
left_join(ci_90.VacUli, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.VacUli <- model_plot_marg_function(coeff.shrub_gradient.VacUli, title_string = "Vaccinium uliginosum", plot_width = 11))
As the quadratic terms for temperature variability and precipitation were not significant, they were removed from the model before re-running:
# new model object with terms removed
write("
model{
# priors
intercept ~ dnorm(0, 0.0001)
b.compet ~ dnorm(0, 0.0001)
b.sri ~ dnorm(0, 0.0001)
b.tri ~ dnorm(0, 0.0001)
b.twi ~ dnorm(0, 0.0001)
sigma.plotgroup ~ dunif(0,100)
tau.plotgroup <- 1/(sigma.plotgroup * sigma.plotgroup)
b.tempjja.x ~ dnorm(0, 0.001)
b.tempjja.x2 ~ dnorm(0, 0.001)
b.tempcont.x ~ dnorm(0, 0.001)
# b.tempcont.x2 ~ dnorm(0, 0.001)
b.precipjja.x ~ dnorm(0, 0.001)
# b.precipjja.x2 ~ dnorm(0, 0.001)
phi ~ dgamma(0.1, 0.1)
# LIKELIHOOD for discrete part
for (i in 1:N_discrete){
cov.dis[i] ~ dbern(mu[i])
logit(mu[i]) <- b_plotgroup[plotgroup.dis[i]] + # ~= random effect of plot group
b.compet * compet.dis[i] +
b.twi * twi.dis[i] +
b.sri * sri.dis[i] +
b.tri * tri.dis[i]
}
# LIKELIHOOD for continuous part
for (j in 1:N_cont){
cov.cont[j] ~ dbeta(p[j], q[j])
p[j] <- mu2[j] * phi
q[j] <- (1 - mu2[j]) * phi
logit(mu2[j]) <- b_plotgroup[plotgroup.cont[j]] + # ~= random effect of plot group
b.compet * compet.cont[j] +
b.twi * twi.cont[j] +
b.sri * sri.cont[j] +
b.tri * tri.cont[j]
}
for (k in 1:N_plotgroups){ # length of total plotgroups
b_plotgroup[k] ~ dnorm(mu.plotgroup[k],tau.plotgroup)
mu.plotgroup[k] <- intercept +
# plot group level predictors, linear and quadratic term
b.tempjja.x * tempjja.tot[k] +
b.tempjja.x2 * (tempjja.tot[k]^2) +
b.tempcont.x * tempcont.tot[k] +
# b.tempcont.x2 * (tempcont.tot[k]^2) +
b.precipjja.x * precipjja.tot[k] # +
# b.precipjja.x2 * (precipjja.tot[k]^2)
}
# add predicted values (derived parameters)
for (m in 1:Nxhat){
phat_compet[m] <- intercept + b.compet * xhat_compet[m]
phat_sri[m] <- intercept +
b.sri * xhat_sri[m]
phat_tri[m] <- intercept +
b.tri * xhat_tri[m]
phat_twi[m] <- intercept +
b.twi * xhat_twi[m]
phat_tempjja[m] <- intercept +
b.tempjja.x * xhat_tempjja[m] +
b.tempjja.x2 * (xhat_tempjja[m]^2)
phat_tempcont[m] <- intercept +
b.tempcont.x * xhat_tempcont[m] # +
# b.tempcont.x2 * (xhat_tempcont[m]^2)
phat_precipjja[m] <- intercept +
b.precipjja.x * xhat_precipjja[m] # +
# b.precipjja.x2 * (xhat_precipjja[m]^2)
for (p in 1:Nxhat2){
phat_tempcontXtemp[m,p] <- intercept +
b.tempcont.x * xhat_tempcont[m] +
# b.tempcont.x2 * (xhat_tempcont[m]^2) +
b.tempjja.x * xhat_tempjja2[p] +
b.tempjja.x2 * (xhat_tempjja2[p]^2)
}
}
}
", file.path("..", "models", "shrub_gradient.VacUli2.jags"))
# specify new set of parameters to be monitored
params_VacUli2 <- c("intercept",
"b.tempjja.x", "b.tempjja.x2",
"b.tempcont.x", # "b.tempcont.x2",
"b.precipjja.x", # "b.precipjja.x2",
"b.compet",
"b.sri",
"b.tri",
"b.twi",
"b_plotgroup[1]","b_plotgroup[2]","b_plotgroup[3]","b_plotgroup[63]",
"sigma.plotgroup",
"phi",
"phat_compet", "phat_sri", "phat_tri", "phat_twi", "phat_tempjja", "phat_tempcont", "phat_precipjja", "phat_tempcontXtemp")
model_out.shrub_gradient.VacUli2 <- jags(shrub_gradient_jags.VacUli.data, # input data
inits = NULL, # JAGS to create initial values
params_VacUli2, # parameters to be saved
model.file = file.path("..", "models", "shrub_gradient.VacUli2.jags"),
n.chains = 3, # no. Markov chains
n.iter = 100000, n.burnin = 70000, # no. iterations & burn-in fraction per chain
n.thin = 2, # thinning rate
DIC = FALSE, # do not compute deviance, pD, and DIC
working.directory = NULL,
progress.bar = "text")
Compiling model graph
Resolving undeclared variables
Allocating nodes
Graph information:
Observed stochastic nodes: 414
Unobserved stochastic nodes: 80
Total graph size: 7819
Initializing model
|
| | 0%
|
|+ | 3%
|
|+++ | 6%
|
|++++ | 9%
|
|++++++ | 11%
|
|+++++++ | 14%
|
|+++++++++ | 17%
|
|++++++++++ | 20%
|
|+++++++++++ | 23%
|
|+++++++++++++ | 26%
|
|++++++++++++++ | 29%
|
|++++++++++++++++ | 31%
|
|+++++++++++++++++ | 34%
|
|+++++++++++++++++++ | 37%
|
|++++++++++++++++++++ | 40%
|
|+++++++++++++++++++++ | 43%
|
|+++++++++++++++++++++++ | 46%
|
|++++++++++++++++++++++++ | 49%
|
|++++++++++++++++++++++++++ | 51%
|
|+++++++++++++++++++++++++++ | 54%
|
|+++++++++++++++++++++++++++++ | 57%
|
|++++++++++++++++++++++++++++++ | 60%
|
|+++++++++++++++++++++++++++++++ | 63%
|
|+++++++++++++++++++++++++++++++++ | 66%
|
|++++++++++++++++++++++++++++++++++ | 69%
|
|++++++++++++++++++++++++++++++++++++ | 71%
|
|+++++++++++++++++++++++++++++++++++++ | 74%
|
|+++++++++++++++++++++++++++++++++++++++ | 77%
|
|++++++++++++++++++++++++++++++++++++++++ | 80%
|
|+++++++++++++++++++++++++++++++++++++++++ | 83%
|
|+++++++++++++++++++++++++++++++++++++++++++ | 86%
|
|++++++++++++++++++++++++++++++++++++++++++++ | 89%
|
|++++++++++++++++++++++++++++++++++++++++++++++ | 91%
|
|+++++++++++++++++++++++++++++++++++++++++++++++ | 94%
|
|+++++++++++++++++++++++++++++++++++++++++++++++++ | 97%
|
|++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
|
| | 0%
|
|*** | 7%
|
|******* | 13%
|
|********** | 20%
|
|************* | 27%
|
|***************** | 33%
|
|******************** | 40%
|
|*********************** | 47%
|
|*************************** | 53%
|
|****************************** | 60%
|
|********************************* | 67%
|
|************************************* | 73%
|
|**************************************** | 80%
|
|******************************************* | 87%
|
|*********************************************** | 93%
|
|**************************************************| 100%
# plot(model_out.shrub_gradient.VacUli2) #check convergence, etc.
Extract coefficients and plot effect sizes:
# extract coefficients
coeff.shrub_gradient.VacUli2 <- model_out.shrub_gradient.VacUli2$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# add 90% CIs
ci_90.VacUli2 <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.VacUli$BUGSoutput$sims.list)-4)){
ci_90.VacUli2[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.VacUli2$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.VacUli2[param, 3] <- names(data.frame(model_out.shrub_gradient.VacUli2$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.VacUli2 <- coeff.shrub_gradient.VacUli2 %>%
left_join(ci_90.VacUli2, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.VacUli2 <- model_plot_marg_function(coeff.shrub_gradient.VacUli2, title_string = "Vaccinium uliginosum", plot_width = 8.5))
# 3x2 grid (vertical layout)
(nuuk_effectsize_plot_grid_ver <- plot_grid(effect_size_plot.BetNan,
effect_size_plot.EmpNig,
effect_size_plot.RhoGro,
effect_size_plot.SalGla,
effect_size_plot.VacUli,
labels = c("a)", "b)", "c)", "d)", "e)"),
label_size = 20,
label_fontface = "plain",
hjust = -2,
ncol = 2))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_betabinom_effect_size_panels_vert.eps"),
# nuuk_effectsize_plot_grid_ver, base_height = 15, base_aspect_ratio = 0.8)
# 3x2 grid reduced predictor terms
(nuuk_reduced_effectsize_plot_grid_ver <- plot_grid(effect_size_plot.BetNan2,
effect_size_plot.EmpNig2,
effect_size_plot.RhoGro2,
effect_size_plot.SalGla2,
effect_size_plot.VacUli2,
labels = c("a)", "b)", "c)", "d)", "e)"),
label_size = 20,
label_fontface = "plain",
hjust = -2,
ncol = 2))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_betabinom_reduced_effect_size_panels_vert.eps"),
# nuuk_reduced_effectsize_plot_grid_ver, base_height = 15, base_aspect_ratio = 0.8)
# 2x3 grid (horizontal layout)
(nuuk_effectsize_plot_grid_hor <- plot_grid(effect_size_plot.BetNan,
effect_size_plot.EmpNig,
effect_size_plot.RhoGro,
effect_size_plot.SalGla,
effect_size_plot.VacUli,
labels = c("a)", "b)", "c)", "d)", "e)"),
label_size = 20,
label_fontface = "plain",
hjust = -2,
ncol = 3))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_betabinom_effect_size_panels_hor.eps"),
# nuuk_effectsize_plot_grid_hor, base_height = 10, base_aspect_ratio = 1.5)
# 2x3 grid reduced predictor terms
(nuuk_reduced_effectsize_plot_grid_hor <- plot_grid(effect_size_plot.BetNan2,
effect_size_plot.EmpNig2,
effect_size_plot.RhoGro2,
effect_size_plot.SalGla2,
effect_size_plot.VacUli2,
labels = c("a)", "b)", "c)", "d)", "e)"),
label_size = 20,
label_fontface = "plain",
hjust = -2,
ncol = 3))
# save plot
# save_plot(file.path("..", "figures", "nuuk_shrub_drivers_betabinom_reduced_effect_size_panels_hor.eps"),
# nuuk_reduced_effectsize_plot_grid_hor, base_height = 10, base_aspect_ratio = 1.5)
unique.compet <- seq(min(SalGla.tot$competC), max(SalGla.tot$competC), by = 0.02)
preds.compet.SalGla <- array(dim=c(length(unique.compet),
nrow(model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list$b_plotgroup)))
for (i in 1:length(unique.compet)){
preds.compet.SalGla[i,] <- plogis(model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list$b.compet[,1]*unique.compet[i])
}
# # preds.continent1biome2 <- array(dim=c(length(unique.compet),
# # length(model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list$b1meanForest[,1])))
# #
# # for (i in 1:length(unique.compet)){
# # preds.continent1biome2[i,] <- plogis(model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list$b1meanSavanna + model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list$brain.x[,2]*unique.compet[i] + model_out.shrub_gradient.SalGla2$BUGSoutput$sims.list$brain.x2[,2]*(unique.compet[i]^2))}
# # continentpreds <- array(c(preds.continent1biome1,preds.continent1biome2), dim=c(length(unique.compet), ncol(preds.continent1biome1),2))
preds.compet.SalGla[1:10,1:10] # [compet,sims]
quantiles.compet.SalGla <- array(NA, dim=c(length(unique.compet),5))
for (i in 1:length(unique.compet)){#compet
#i<-1
quantiles.compet.SalGla[i,]<-quantile(preds.compet.SalGla[i,], c(0.025, 0.25, 0.5, 0.95, 0.975))
}
quantiles.compet.SalGla[1:5,1:5]
predsmatrix.compet.SalGla <- data.frame(expand.grid(competC=c(unique.compet)))
predsmatrix.compet.SalGla$l95CI<-rep(NA, nrow(predsmatrix.compet.SalGla))
predsmatrix.compet.SalGla$l90CI<-rep(NA, nrow(predsmatrix.compet.SalGla))
predsmatrix.compet.SalGla$med<-rep(NA, nrow(predsmatrix.compet.SalGla))
predsmatrix.compet.SalGla$u90CI<-rep(NA, nrow(predsmatrix.compet.SalGla))
predsmatrix.compet.SalGla$u95CI<-rep(NA, nrow(predsmatrix.compet.SalGla))
head(predsmatrix.compet.SalGla)
# # rainfallseq<-rep(1:length(unique.compet),each=2)
for (i in 1:nrow(predsmatrix.compet.SalGla)){
predsmatrix.compet.SalGla[i,c(2:6)] <- quantiles.compet.SalGla[i,]
}
predsmatrix.compet.SalGla$compet <- predsmatrix.compet.SalGla$competC + attr(scale(SalGla.tot$compet), 'scaled:center')
(pred_plot.compet.SalGla <- ggplot() +
# competition is modeled on plot level, so use input data as is
geom_point(data = SalGla.tot,
aes(x = compet,
y = cover),
size = 2,
position = position_jitter(width=0, height=.01),
alpha = 0.5) +
# plot line of predicted median values
geom_line(data = predsmatrix.compet.SalGla,
aes(x = compet,
y = med),
colour = "orange",
alpha = 1,
size = 3) +
# scale_colour_manual(values=c("darkgreen","orange","darkgreen","orange"),name="Biome") +
# plot predicted 95% CI
geom_ribbon(data = predsmatrix.compet.SalGla,
aes(x = compet,
ymin = l95CI,
ymax = u95CI),
fill = "orange",
alpha = 0.2) +
# plot predicted 90% CI
geom_ribbon(data = predsmatrix.compet.SalGla,
aes(x = compet,
ymin = l90CI,
ymax = u90CI),
fill = "orange",
alpha = 0.4) +
# scale_fill_manual(values=c("darkgreen","orange"),name="Biome") +
labs(x = "summed abundance of taller-growing species per plot",
y = "rel. no. hits per plot") +
theme_bw())
looks weird; prediction does not relate to data points in x or y range.
Let’s run another model with some derived values to see if we get the same numbers out as the way we have extracted them above (i.e., by scaling back). That will tell us whether it’s a problem with the back-scaling or a more fundamental issue.
# assemble data: SalGla as example ----
shrub_gradient_jags.SalGla.xhat.data <- list(
# plot level predictors, for discrete...
cov.dis = SalGla.dis$cover,
plotgroup.dis = SalGla.dis$plotgroup.NUM, #AB added this
# isocline.dis = SalGla.dis$site_alt.NUM,
# inclin_down.dis = SalGla.dis$inclin_downC,
sri.dis = SalGla.dis$sriC,
tri.dis = SalGla.dis$triC,
twi.dis = SalGla.dis$twiC,
compet.dis = SalGla.dis$competC,
N_discrete = nrow(SalGla.dis),
# ...and continuous part of the data
cov.cont = SalGla.cont$cover,
plotgroup.cont = SalGla.cont$plotgroup.NUM, #AB added this
# isocline.cont = SalGla.cont$site_alt.NUM,
# inclin_down.cont = SalGla.cont$inclin_downC,
sri.cont = SalGla.cont$sriC,
tri.cont = SalGla.cont$triC,
twi.cont = SalGla.cont$twiC,
compet.cont = SalGla.cont$competC,
N_cont = nrow(SalGla.cont),
# plot group level predictors
tempjja.tot = SalGla.tot$tempjja_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)], # one value per tXpg
# tempmax.tot = SalGla.tot$tempmax_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
# tempmin.tot = SalGla.tot$tempmin_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
tempcont.tot = SalGla.tot$tempcont_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
precipjja.tot = SalGla.tot$precipjja_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)],
# precipjfmam.tot = SalGla.tot$precipjfmam_ts_30C[!duplicated(SalGla.tot$plotgroup.NUM)]
N_plotgroups = length(unique(SalGla.tot$site_alt_plotgroup_id)),
xhat = seq(from = min(SalGla.tot$tempcont_ts_30C), to = max(SalGla.tot$tempcont_ts_30C), length.out = 100),
Nxhat = 100
# # site/alt level predictors
# alt.tot = SalGla.tot$altC[!duplicated(SalGla.tot$site_alt.NUM)],
# N_isoclines = length(unique(SalGla.tot$site_alt_id))
)
str(shrub_gradient_jags.SalGla.xhat.data)
We’ll add some predicted values (phat) to the model:
Specify the parameters to be monitored:
Run model:
Extract coefficients and make graph of data and prediction curve:
# extract coefficients
coeff.shrub_gradient.SalGla2.xhat <- model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.SalGla2.xhat),"[[]",fixed=FALSE), "[", 1))) #%>% print
# # add 90% CIs
# ci_90.SalGla2.xhat <- data.frame(q5 = NA,
# q95 = NA,
# param = NA)
# for (param in 1:(length(model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$sims.list))){
# ci_90.SalGla2.xhat[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$sims.list[param])[,1],
# probs = c(0.05, 0.95))
# ci_90.SalGla2.xhat[param, 3] <- names(data.frame(model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$sims.list))[param]
# }
# join to coefficients table
coeff.shrub_gradient.SalGla2.xhat <- coeff.shrub_gradient.SalGla2.xhat %>%
# left_join(ci_90.SalGla2.xhat, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
# l90 = q5,
# u90 = q95,
u95 = "97.5%",
Rhat) # %>% print
# assemble predicted and predictor values
phats <- coeff.shrub_gradient.SalGla2.xhat %>%
# filter for predicted values
filter(param %in% c(paste0("phat[", seq(from = 1, to = 100), "]"))) %>%
# add xhats column
mutate(xhats = seq(from = min(SalGla.tot$tempcont_ts_30C),
to = max(SalGla.tot$tempcont_ts_30C),
length.out = 100))
# back-center and back-scale predictor values (xhats)
phats$tempcont <- phats$xhats*attr(scale(SalGla.tot$tempcont_ts_30), 'scaled:scale') + attr(scale(SalGla.tot$tempcont_ts_30), 'scaled:center')
ggplot() +
# tempcont is modelled at plotgroup level, so reduce base data (points layer) to plotgroup level
geom_point(data = SalGla.tot %>% group_by(site_alt_plotgroup_id) %>% summarise(tempcont = mean(tempcont_ts_30), cover = mean(cover)),
aes(x = tempcont,
y = cover),
size = 2,
position = position_jitter(width=0, height=.01),
alpha=0.5) +
# draw line of predicted values
geom_line(data = phats,
aes(x = tempcont,
y = plogis(mean)),
colour = "orange",
alpha = 1,
size = 3) +
# scale_colour_manual(values=c("darkgreen","orange","darkgreen","orange"),name="Biome") +
# draw predicted 95% CI
geom_ribbon(data = phats,
aes(x = tempcont,
ymin = plogis(l95),
ymax = plogis(u95)),
fill = "orange",
alpha = 0.2) +
# draw predicted 90% CI
# geom_ribbon(data = phats,
# aes(x = tempcont,
# ymin = plogis(l90),
# ymax = plogis(u90)),
# fill = "orange",
# alpha = 0.4) +
# scale_fill_manual(values=c("darkgreen","orange"),name="Biome") +
labs(x = "annual temperature variability [°C]",
y = "rel. no. hits per plot") +
theme_bw()
Extract coefficients and plot effect sizes:
# load model object
load(file = file.path("..", "data", "processed", "model_output_SalGla2_with_derived_values.RData"))
# extract coefficients
coeff.shrub_gradient.SalGla2.xhat <- model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$summary %>%
as.data.frame %>%
select('mean','sd','2.5%','97.5%','Rhat') %>%
# add identifying info to data frame
rownames_to_column(var = "param")
# mutate(param = as.vector(sapply(strsplit(rownames(coeff.shrub_gradient.SalGla),"[[]",fixed=FALSE), "[", 1))) #%>% print
# add 90% CIs
ci_90.SalGla2.xhat <- data.frame(q5 = NA, q95 = NA, param = NA)
for (param in 1:(length(model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$sims.list)-4)){
ci_90.SalGla2.xhat[param,1:2] <- quantile(data.frame(model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$sims.list[param])[,1], probs = c(0.05, 0.95))
ci_90.SalGla2.xhat[param, 3] <- names(data.frame(model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$sims.list))[param]
}
# join to coefficients table
coeff.shrub_gradient.SalGla2.xhat <- coeff.shrub_gradient.SalGla2.xhat %>%
left_join(ci_90.SalGla2.xhat, by = "param") %>%
# reorder and rename cols
select(param, mean, sd,
l95 = "2.5%",
l90 = q5,
u90 = q95,
u95 = "97.5%",
Rhat) %>% print
# effect size plot
(effect_size_plot.SalGla2.xhat <- model_plot_sig_function(coeff.shrub_gradient.SalGla2.xhat, title_string = "Salix glauca", plot_width = 11))
extract predicted values:
# get statistical summary
model_xhat.out.df <- as.data.frame(model_out.shrub_gradient.SalGla2.xhat$BUGSoutput$summary[,c('mean','sd','2.5%','97.5%','n.eff','Rhat')])
model_xhat.out.df$Param <- sapply(strsplit(rownames(model_xhat.out.df), "[[]", fixed=F),"[",1)
p.out <- model_xhat.out.df[model_xhat.out.df$Param=="phat",]
p.out$xhat <- shrub_gradient_jags.SalGla.xhat.data$xhat #add the "x" values you used to predict phats
plot along with back-scaled values:
pred_plot.compet.SalGla +
geom_line(data = p.out,
aes(x = xhat,
y = mean),
colour = "darkgreen",
alpha = 1,
size = 2) +
geom_ribbon(data = p.out,
aes(x = xhat,
ymin = `2.5%`,
ymax = `97.5%`),
alpha = 0.3,
fill = "green")
# save.image(file = file.path("..", "data", "processed", "nuuk_shrub_drivers_output_nonsigx2removed.RData"))
# load ws for further analyses/plotting etc)
# load(file.path("data", "processed", "nuuk_shrub_drivers_output_nonsigx2removed.RData"))